Deriving a New Crustal Model of Northern Adria: The Northern Adria Crust (NAC) Model

We presented a new 3D model of the geophysical properties of the crust (namely depth of the Moho and VP, VS, density, Young’s modulus, and shear modulus) of the northern tip of the Adria microplate that we called NAC (Northern Adria Crust). The horizontal dimensions of the physical properties variations are optimized at 5 × 5 km and the vertical dimension at 1 km. NAC has been built by critically choosing and integrating all available information about the depth of the main interfaces and the physical properties of the crust. We started from a VP dataset, and we converted it in VS and density by using empirical relations, tuned through the comparison with the available data from local tomographic inversion, and taking into account the lithologies of the area. Uncertainties and reliability of the model are quantified, taking into account the data quality and the interpolation procedure. NAC has two versions, different in the structure of the Moho interface: the first considers one continuous surface for the whole area, while the second implies three separate surfaces for the Adria microplate, Eurasia, and the Pannonian fragment. The differences between the two models are minimal, but the available data better sustain the solution of the fragmented crust. For its characteristics of multiparametric information and resolution, NAC can be precious for any purpose and use where a detailed knowledge of the crustal structure of this area is required. Moreover, it is easy to improve NAC, including new information on the crustal structures, when they will be available.


INTRODUCTION
An accurate model of the crust is the basis for modeling in various fields as seismology, geodynamics, gravity (e.g., Zhu and Tromp, 2013;Faccenna et al., 2014;Métois et al., 2015;Root et al., 2015). Highly detailed modeling is of particular importance in areas of continental collision, i.e., the orogenic belts, still offering unanswered questions on the geometry of the boundaries, the reciprocal roles, and the relative motions of the converging plates. One of these areas is the area across North-Eastern Italy, Austria, Western Slovenia, and Croatia: namely, the transition between Eastern Alps and External Dinarides. It occupies the northernmost edge of the convergent margin between Eurasia and the Adria microplate (Figure 1), which is relatively aseismic, and encircled by active orogenic belts (Argand, 1924;McKenzie, 1972;Channell et al., 1979; FIGURE 1 | Map that shows the area of the model (red rectangle). The origin point of the reference system of the model (O) is located at 10.2E, 45.6N, and the model extends 380 km to the East, 300 km to the North, and from 4 km a.sl. to 60 km b.s.l in depth. TW -Tauern Window. Blue line: main tectonic lineaments: PAL -Periadriatic Lineament; GL -Giudicarie line. Black lines: boundaries between Moho fragments, according to Brückl et al. (2007). Anderson and Jackson, 1987). Seismological and geodetic studies helped to establish that the Adria microplate counterclockwise rotates relatively to Eurasia around a pole of rotation located in the western Po Plain (e.g., Battaglia et al., 2003;D'Agostino et al., 2005D'Agostino et al., , 2008Serpelloni et al., 2005).
The seismic investigations in the Eastern and Southern Alps started in the early sixties, and developed in the seventies of the last century, with Deep Seismic Soundings (DSS) profiles, later revised by Scarascia and Cassinis (1997). In the last years, new information about the crustal properties of the area became available thanks to a series of controlled-source seismic experiments: TRANSALP (TRANSALP Working Group, 2002), ALP 2002 Grad et al., 2009a;Šumanovac et al., 2009), and CROP (Finetti, 2005) experiments.
One of the main reasons for such high interest is the vivid debate about the subduction directions and the necessity of discriminating between different tectonic models. The high-resolution TRANSALP transect in the area of the Tauern Window (about 12 • E) and the complementary experiments provided evidences of the subduction of the Eurasia below the Adria microplate (TRANSALP Working Group, 2002;Kummerow et al., 2004;Lüschen et al., 2004Lüschen et al., , 2006Bleibinhaus and Gebrande, 2006;Castellarin et al., 2006). While confirming the southward subduction to the west, data from teleseismic tomography indicate a change of the subduction direction to the East (Lippitsch et al., 2003;Schmid et al., 2004;Kissling et al., 2006;Handy et al., 2015). A different perspective was offered by the CELEBRATION 2000, ALP 2002 WAR/R, and ALPASS data in the transition area between the Alps, Dinarides and Carpathians, interpreted in terms of a Moho fragmentation between Eurasia, the Adria microplate, and the Pannonian fragment (Behm et al., 2007;Brückl et al., 2007Brückl et al., , 2010Šumanovac et al., 2009;Mitterbauer et al., 2011;Figure 1), so that the subduction polarity of Eurasia plate would be unique in the whole Alpine region. To the south of the triple junction between the three fragments, the Adria lithosphere would underthrust below the Dinarides and the Pannonian fragment (Brückl et al., 2010;Šumanovac et al., 2016).
Recently, also ambient noise tomography was applied in the studied region (Molinari et al., 2015b;Behm et al., 2016;Guidarelli et al., 2017). Kästle et al. (2018), through the joint inversion of ambient noise and earthquake phase velocity measurements, confirmed the differences in the velocity structure of the Central and Eastern Alps, but could not give a definitive answer to the questions about the subduction direction. Kästle et al. (2020) reviewed the different slab break-off models, offering a new combined interpretation of body-wave and surface-wave tomography, to be verified by the forthcoming data. Moreover, two of the complementary experiments of the AlpArray program -EASI (AlpArray Seismic Network, 2014;Hetényi et al., 2018) and SWATH-D (Heit et al., 2017) -focus on this area. The two seismological profiles spaced 15 km apart of EASI revealed the complexity of the crustal structure (both at shallow and deep levels) while supporting the model of an Adria subduction below Eurasia, with a steep northwards dipping slab (Hetényi et al., 2018). New insights on the crustal structure, as, e.g., done by Spooner et al. (2019), are expected to add on in the next years from the SWATH-D data analysis and other AlpArray complementary experiments and continuation.
We focus on the region to the south of the PeriAdriatic Lineament (PAL in Figure 1). Here, as a result of the counterclockwise rotation of the Adria microplate, the convergence rate increases across the Alps from W to E, with a maximum N-S shortening of about 2 mm/a over ∼80 km at the easternmost part of the Southern Alps, corresponding to the Italian Friuli seismic region (D'Agostino et al., 2005;Grenerczy et al., 2005;Völksen et al., 2018). The latter is the seismically most active area of the Alps, with several destructive earthquakes in historic times, including the M = 6.4 1976 Friuli earthquake (e.g., Burrato et al., 2008;Santulin et al., 2018, and references therein). A local seismometric network monitors the region since 1977, complemented, since 2002, by a GNSS network (Bragato et al., 2013;Zuliani et al., 2018), making the region a sort of natural laboratory for the study of complex tectonic processes. The seismicity is mainly located in the upper crust (Viganò et al., 2015;. The stress and strain tensor inversions from the focal mechanisms reveal in the region a complex stress field with a resulting heterogeneous deformation pattern, with variations on a scale of tens of km (Bressan et al., 2018a). The inversion of fault plane solutions reveals a compression ranging from west to east from NW-SE to N-S, and NNE-SSW, in agreement with the counter-clockwise motion and interaction of the Adria microplate with the Eurasia plate, and the geodynamic frame described above (Bressan et al., 2018a). The spatial seismicity pattern is conditioned, especially in the shallowest layers (0-10 km), by the crustal heterogeneities .
Recently, the increasing number and quality of GNSS data provided new insight into the kinematics of the area (Caporali et al., 2009;Métois et al., 2015;Rossi et al., 2016;Serpelloni et al., 2016). The observations confirmed the shortening of about 2 mm/a across the southern front of the Eastern Alps, and the increasing eastward rotation of velocities toward the Pannonian Basin, whereas an N-S directed extension can be related to a gravitational collapse within the Tauern Window. The sharp boundaries between the strain velocities observed are in agreement with the triple junction between Eurasia, Adria microplate, and the Pannonian fragment (Möller et al., 2011;Völksen et al., 2018).
A finite element modeling can help to relate the seismic activity with the crustal deformation observed by the GNSS measurements through the crustal geophysical properties and analyze how the forces at plate boundaries influence the intraplate stresses. The 2D finite element model presented by Bada et al. (1998), although mainly focused on the Pannonian basin, allowed to explain most of the stress patterns observed by modeling the forces acting at the plate boundaries. Rossi et al. (2005) 3D finite element model of NE Italy related the seismic focal mechanisms and the seismicity distribution to the counterclockwise rotation of the Adria microplate relative to Eurasia.
Our goal is constructing a more complex 3D model of the northern tip of the Adria microplate, based on the actual knowledge about the crust structure, and physical properties, enabling modeling also the stresses induced by density contrasts in the crust (Ranalli, 1992;Heidbach et al., 2007), horizontal gradients of gravitational potential energy (Marotta and Splendore, 2014;Métois et al., 2015), topography (Heidbach et al., 2007), and plate boundaries (Richardson, 1992) in this critical area of the Alps. Detailed knowledge of the elastic moduli (derivable from the seismic velocities and density) distribution in the crust is, hence, necessary. Moreover, mapping in detail the Moho discontinuity can help to verify its possible fragmentation while providing keys for understanding the dynamic processes causing crustal growth, accretion, delamination, and underplating (Carbonell et al., 2013, and references therein).
The physical properties (namely: seismic velocities, rock bulk density, and elastic moduli) of the upper crust of the Friuli seismic region were studied through laboratory measurements on rock samples, representative of the most common lithologies, sonic logs, local earthquake tomography, and seismo-gravimetric inversion (Faccenda et al., 2007;. Local earthquake tomography was also applied in other parts of the Southern Alps (Anselmi et al., 2011;Viganò et al., 2013).
The crustal models at the continental scale developed in the last 10 years include in whole or in part, the most recent data about the crustal structure. These datasets include (from old to young) EuCrust-07 (Tesauro et al., 2008), the Moho map of Grad et al. (2009b) (ESC Moho in the following), EPCrust (Molinari and Morelli, 2011), the Moho maps of Italy by Di Stefano et al. (2011) and Spada et al. (2013), and EUNAseis (Artemieva and Thybo, 2013). A density model of the Alpine lithosphere, constrained by 3D gravity modeling and based on all the recent seismic experiments is provided by Spooner et al. (2019).
EuCrust-07, EPCrust, and EUNAseis models share some critical features like a uniform parameterization that allows a simple conversion in a numerical mesh. They are an essential reference for constructing a regional model, but, because of their scale and resolution, they cannot be directly used for studies focused on the analysis of small scale features (tens of km). More detailed models, developed for seismological studies, cover only marginally the region of our interest (Vuan et al., 2011;Molinari et al., 2012Molinari et al., , 2015a. The present study fills this partial gap, by presenting a new 3D crustal model of the northern tip of the Adria microplate, that we call NAC (Northern Adria Crust). It has been constructed by critically collecting and integrating all available information in a 3D model (from the Moho to the topographic surface) of the crustal properties as illustrated in section Data. In section Model Construction we explain the procedure adopted to integrate the various data, interpolating them on the model grid, and evaluating the uncertainty of the interface depth and the physical properties (seismic velocities, density, Young, and shear modulus) that we retrieve. The uncertainty of the interfaces' depth and on the various physical properties is quantified by taking into account the data quality and coverage and the effects of the interpolation procedure. In section Physical Properties of the Crust, we critically analyze the fit of different empirical relationships to convert P-wave velocity (V P ) into density (ρ) and S-wave velocity (V S ). We compared the results with the available information from independent measurements of V S , laboratory experiments, and borehole measurements. The uncertainty estimation, as well as the considerations on the physical properties empirical relationships can be of interest for applications in other parts of the world and tectonic frameworks. The resulting 3D crustal properties multiparameter model is described in section Main Features of the Model, and discussed in section Discussion. It is defined on a 5 × 5 km grid (∼0.05 • × • 0.05), enabling the modeling of the observed stress/strain variations. It can be, therefore, a useful tool for future local scale studies in the fields of gravity, geodesy, geodynamics, and seismology. Moreover, NAC's structure is such that it will be easy to include new information and update the model.

DATA
The area here considered includes the eastern Southern Alps (to the east of the Giudicarie line and the south of the PeriAdriatic Lineament -GL and PAL, respectively in Figure 1), the southern sector of the Eastern Alps, the north-western external Dinarides, the Venetian and Friulian plains, part of the Po plain, and the northern part of the Adriatic Sea (Figure 1). We collected information about the geophysical properties (seismic velocities and bulk rock density) of the crustal layers and their geometry (Figure 2). In this section, we describe the data sources that we used in the construction of the model ( Table 1).
The criterion adopted to choose the data to be included in the model, as well as the way of doing it, is the data uncertainty. Uncertainties of the Moho map of Behm et al. (2007) in the area of our interest are around 2 km (Behm, 2006). Brückl et al. (2007), Grad et al. (2009a), andŠumanovac et al. (2009) reported similar uncertainties of the interfaces depth: 2-3 km for mid-crustal boundaries and 1-2 km for the Moho. For the other seismic profiles, we followed the estimations of Grad et al. (2009b) that assigned an uncertainty of 6-8% to the Moho depth retrieved from old seismic profiles (i.e., about 3 km for a Moho depth of 40 km). To have a conservative estimate of the uncertainties, we used the upper bound of these ranges. We compared the various seismic profiles at the intersections (Supplementary Figure S1), and we found that the average difference between the Moho depth is about 3.6 km. The most significant difference (about 10 km) is observed between the data of the CEL10/Alp04 and the profile ALP75 , where Brückl et al. (2010) trace the boundary between Eurasia and the Pannonian fragment.
The teleseismic receiver functions have a lower resolution than the seismic profiles. Therefore, we excluded from our dataset the stations from Orešković et al. (2011);Bianchi et al. (2015), and Šumanovac et al. (2016) that cover the same areas as the seismic profiles, retaining the ones that provide information about the Moho depth in zones less covered by seismic profiles, as the South-Western part of the model (Figure 3). We used only the stations of Piana Agostinetti and Amato (2009), and of Miller and Piana-Agostinetti (2012) of quality classes 1 and 2. We also excluded the Moho depth data with a relative uncertainty >20%. When for the same station, more than one estimation of Moho depth is available, we checked if they are compatible and took the one with the lowest error. We did not include data from Hetényi et al. (2018) since they did not report a unique crust-mantle boundary, i.e., a clear Moho, and the proposed depth shows significant differences with the seismic profiles (Alp01, Alp02, and three profiles from Scarascia and Cassinis, 1997) and with other receiver function results in the same area (Bianchi et al., 2015).
Seismic profiles Bleibinhaus and Gebrande, 2006;Behm et al., 2007;Brückl et al., 2007;Grad et al., 2009a;Šumanovac et al., 2009) are the principal sources of information also about seismic velocities (mainly V P ). We followed Artemieva and Thybo (2013) to assign V P accuracy to the different seismic profiles. Local earthquake tomography studies provide information about V P and V S in the Southern Alps and External Dinarides (Anselmi et al., 2011;Viganò et al., 2013). Gravimetric studies at regional scale provide information on the density of the crustal layers (Braitenberg et al., 1997;Cassinis et al., 1997;Dal Moro et al., 1998;Ebbing et al., 2001Ebbing et al., , 2006Šumanovac et al., 2009;Šumanovac, 2010). The geophysical properties of the upper crust in the area between the Southern Alps and the External Dinarides were also studied by laboratory measurements (Faccenda et al., 2007) as well as through the inversion of seismic and gravimetric data (de Franco et al., 2004;. As explained more in detail in the following sections, the observations on ρ and V S do not enter directly in the NAC construction but are used to select the best empirical relation between V P and the other geophysical parameters and to check the results. Data about the shallowest part of the crust are completed from other sources. A structural model of the Po and Venetian plain of Vuan et al. (2011) provides the depth of the bottom of the sedimentary layer. We also considered the stratigraphy from boreholes Cimolino et al., 2010;ViDEPI, 2015) in the Venetian and Friulian plain, and shallow seismic profiles (Fantoni et al., 2003;Nicolich et al., 2004). We defined the uncertainty of the depth derived from the model of Vuan et al. (2011) by comparing it with independent inputs (boreholes). We included the outcropping of sedimentary rocks at the border of the plains and on the eastern side of the Adriatic Sea. We also added the thickness of the soft sediments in the Klagenfurt (Nemes et al., 1997) and Ljubljana basins (Gosar and Lenart, 2010;. The geophysical properties of the sediments are supplied by the seismo-gravimetric inversion of Tondi et al. (2019) for the Po Plain, the seismological model of Vuan et al. (2011) for the Venetian plain and the Po plain and by Giustiniani et al. (2008) for the Friulian plain.

MODEL CONSTRUCTION
NAC covers an area that extends from the origin point (at 10.2 • E and 44.6 • N) for 380 km to the east and 300 km to the north, hence, covers the region 10.
The model construction is meant to include new information quickly and to convert it in a mesh for numerical computations conveniently. To this aim, a uniform parameterization is desirable, and the procedure used in the model construction should be automatized as much as possible. The choice of the parameters of NAC, described in the following, is similar to what Kelly et al. (2007) did for the UK and Ireland. Our model comprises a first layer of sediments (S), and the crust (C), and is defined by three interfaces: the topography (T), the bottom of the sedimentary layer (BS), and the Moho. Due to the remarkable difference in depth recognized between the Eurasia and the Adria Moho surfaces Kummerow et al., 2004;Brückl et al., 2007;Hetényi et al., 2018) the use of a unique continuous surface representing the Moho could be inappropriate and produce inaccurate results in future modeling. Hence, we consider two versions of the Moho: a unique continuous surface (MO1) and a composite Moho (MO2), consisting of three separated fragments: (i) the Adria microplate (AD); (ii) the Eurasia (EU); (iii) the Pannonian fragment (PA). Compared with other recent crustal models (EuCRUST-07; EPcrust), we adopted a single layer for the crust, in which the seismic velocities (V P , V S ), the density (ρ), and the elastic parameters (Young modulus, E, and shear modulus, µ) are variable, both laterally and vertically.
The physical properties variations are sampled on a grid with cells of 5 × 5 km horizontally, and a vertical dimension of 1 km. The elastic moduli were derived from the values V P , V S , and ρ. In the absence of V S and ρ, we started from a V P dataset, and we converted it into V S and ρ by using empirical relations (section Physical Properties of the Crust). We applied to V S and ρ datasets the same interpolation procedure described in the following for V P data.
The input data is processed consistently in the whole model. The profiles described in the previous section were digitized from existing literature with a regular spacing of 6 km, similar to the estimated average horizontal resolution obtained from controlled-source seismology (CSS) methods for the Moho interface (Waldhauser et al., 1998). We took account of differences at intersections between 2-D profiles by averaging the values of Moho and velocities and, if necessary, by increasing the uncertainties assigned to the data. In the absence of an explicit interpretation of the authors, we took the Moho as the upper boundary where the V P turns to values >7.8 km/s. We extracted the well-resolved Moho depths from the Moho map of Behm et al. (2007) on their regular grid of 20 × 20 km. We defined separated datasets for the three Moho domains to define the MO2. In the eastern part of NAC, we used the plate boundaries defined by Brückl et al. (2010) through their elastic plate modeling. In the central and western parts, we considered the northern edge of the Adria Moho reported by Cassinis (2006) who compared the data from Scarascia and Cassinis (1997) with the results of the TRANSALP profile. The 2D velocity model of each seismic section was converted into a set of 1D velocity profiles. We took the V P values from local earthquake tomography only for those nodes identified by  and Anselmi et al. (2011) as reliable according to a threshold of the spread function (a measure of the blurring of the tomographic image) and by Viganò et al. (2013) according to checkerboard, recovery, and synthetic data tests. The data was transformed from the geographical coordinates to a Cartesian coordinate system to obtain a regular spacing in the area and to convert the model in a mesh suitable for numerical modeling. We used the UTM (Universal Transverse Mercator) projection (zone 33N). The coordinate manipulation was done using pyproj (PROJ contributors, 2018).
As highlighted by Artemieva and Thybo (2013), the resolution of the final model depends mainly upon the input data coverage. Therefore, we analyzed the spatial distribution of the data collected on the depth of each interface and the V P at two different depths (5 and 25 km) by computing the minimum distance (horizontal for the interfaces and considering also the depth for V P ) between the points inside the region (on a regular grid with step of 1 km) and the location of the data source (d min ) (Figure 3). For MO1, 70% of the surface of our region has d min (20 km from the data points, while 90% of the surface of our region has dmin (40 km from them. The area less covered by the data is the South-Western corner. Only a little part of the Adriatic Sea has dmin > 20 km from the data points used for BS. The dmin distribution is are very similar for the two depths and gives a first approximation of the reliability of the final model in the different parts of the model. We interpolated the data points into a regular grid (2D for interfaces and 3D for the velocity model) by kriging using the gstat package (Pebesma, 2004;Gräler et al., 2016), that was used also for the variogram analysis and cross-validation.
Kriging is a collection of generalized linear regression techniques for interpolating sparse geospatial data (Olea, 1999;Oliver and Webster, 2014). The interpolated values are estimated from a weighted sum of sample points using weights that make it unbiased and minimize the variance of the estimation error. Unbiasedness and its optimality in a minimum mean square error sense made kriging a best linear unbiased estimators (Olea, 1999). We introduce some definitions following Olea (1999). The samples are assumed to be a partial realization of a random function Z (x) of the location x. An important assumption of kriging is that the spatial autocorrelation of the variable is known in the form of the semivariogram or covariance: this information is used to estimate the weights of the data points by minimizing the variance (Var) of the estimation error. Given two locations x and x + h inside the field of Z with a constant mean, the semivariogram γ h is: According to this definition, the semivariogram is not a single number but a continuous function of a variable h, called the lag. The lag is a vector; therefore, the semivariogram depends not only on the magnitude of the separation but also on the azimuth of the line through the pairs. An assumption is that the semivariogram is independent of location and depends only on the separation of the pair of locations considered (stationarity). The lag at which the semivariogram reaches a constant value is called the range. The value of the semivariogram beyond the range is named the sill. The values of Z in two points separated by a distance larger than the range are stochastically independent. An estimate of the semivariogram is also called an experimental semivariogram. The semivariogram is considered to be isotropic when variations in the azimuth do not produce significant changes in the experimental semivariogram; otherwise, it is called anisotropic. The experimental semivariogram cannot be used directly, and, therefore, is replaced by a semivariogram model.
As explained in the following, we used both the ordinary kriging (Matheron, 1965) and the universal kriging (Matheron, 1969). The ordinary kriging can be applied if the random function has a constant (but unknown) value. The universal kriging is a generalization that does not require a constant mean. In this case, the random function has a drift: a gentle underlying fluctuation in the regionalized variable. The drift m( x) is the expected value E[:] of the random function: A polynomial drift model is the summation The semivariogram, γ Y , is defined on the residuals of the drift (Matheron, 1969). For both ordinary and universal kriging, the weights are derived from a constrained optimization that is solved through the Lagrange method of multipliers. The estimation at a site − → x 0 is given by: where − → x i are sampled sites and λ i are the kriging weights. The error variance of the estimation for ordinary kriging σ 2 OKI is given by: where µ l is the Lagrange multiplier resulting from the solution of the ordinary kriging system. The error variance of the estimation for ordinary kriging σ 2 UKI is given by: where µ l are the Lagrange multipliers and f l are the terms in the definition of the polynomial drift. For each input dataset (depth of MO1, MO2, BS, and V P field), we proceeded as follows. First, we looked for any evident drift (spatial trend). This is the case for V P along the vertical direction: therefore, we defined a linear model for it. The V P residuals relative to the drift are computed through an ordinary least square estimate of the drift (Equation 2). Then, we computed the experimental semivariograms along different directions to identify the possible anisotropy (Figure 4). The directional semivariograms for the V P residuals relative to the drift have a range in the horizontal plane five times greater than in the vertical direction. Therefore, we included the anisotropy in the semivariogram model for V P . We fitted different models (Gaussian, spherical, and exponential) at each experimental semivariogram, and we chose the best one by cross-validation, comparing the mean error, the mean squared error and the mean squared deviation ratio, which is the mean of the squared errors, divided by the corresponding kriging variances (Oliver and Webster, 2014). Because in all cases we obtained similar results using the exponential and spherical models, whereas the Gaussian model provided the worst results, we used the exponential model for MO1 and each fragment of MO2 and spherical model for V P . The sill obtained for the exponential model of MO1 variogram has a value of ∼80 km, while the one for AD is ∼50 km and the ones for EU and PA are ∼20 km. The short distance variability (nugget) of AD, EU, and PA are similar and lower than the one of MO1. Because the semivariograms for EU and PA fragments are significantly different from that of AD, the stationarity assumption is violated considering a single surface (MO1). The selection of the optimal grid resolution should consider the average spacing between the closest point pairs and the spatial correlation structure (a variable that is spatially auto-correlated at shorter distances would require higher resolution and vice versa) of the input dataset (Hengl, 2006). Since in the horizontal plane the average spacing between the nearest points pairs for MO1 (and MO2) and V P is about 5 km, we defined an XY grid step of 5 km, that is near to the sampling spacing used for seismic profiles. Because the range in the vertical direction is five times smaller than in the horizontal plane, we adopted a step of 1 km along the Z-axis.
To estimate the kriging weights, we also considered the inverse of the square of data uncertainties (Pebesma, 2001). We used the ordinary kriging for the depths of MO1, MO2, and BS, with stationary mean, while the universal kriging for V P data, since it shows a trend with depth. We filtered the interpolated values using a 10 km wide Gaussian filter Smith, 1991, 1998;Wessel et al., 2013) to remove the information not supported by the grid spacing.
To construct MO2, we independently interpolated the datasets for the three segments on three separate surfaces defined by the boundaries described before, and then we combined them in a single surface.
We interpolated the crustal V P values on the regular grid nodes of the model between the Moho and the topographic surface or (where it is present) the bottom of the sedimentary layer.
The uncertainty in the final model (σ TOT ) is composed of three terms: the first derives from the uncertainty in the input data (σ UNC ), the second from the interpolation process (σ INT ) and, only for the physical properties derived from V P , the dispersion of the used empirical relations (σ EMP ). We combined them using the law of error propagation assuming that the uncertainty associated with each term conforms to a Gaussian distribution: We computed the first term interpolating the uncertainty of the input data with the same procedure used for the data itself. Kriging directly provided the second term of uncertainty (σ 2 OKI or σ 2 UKI ). For the physical properties derived from V P , we propagated the uncertainties of original V P values, including also the dispersion of the used empirical relations, computing the partial derivatives of the empirical relations with respect to V P .
As described, the uppermost interface of the NAC model is constituted by the topography. Because the input data for the topography ETOPO1 (Amante and Eakins, 2009) is sampled at 1 arc-minute (about 1.3 km along longitude and 1 along latitude), we resampled them on the coarser grid of NAC, using the adjustable tension continuous-curvature surface gridding algorithm (Smith and Wessel, 1990); the resulting grids were then spatially filtered using a 5 km wide Gaussian filter Smith, 1991, 1998;Wessel et al., 2013).

PHYSICAL PROPERTIES OF THE CRUST
The data about V S in the region is sparse and uneven. Hence, we chose to use empirical functions to convert V P into V S (Brocher, 2005), selected, based on the measured data in the different parts of the model. Brocher (2005) reported two different relations between V P and V S : the "Brocher's regression fit" (BRF) for all the lithologies except calcium-rich and mafic rocks, gabbros, and serpentinites and for 1.5 < V P < 8 km/s and the "mafic line" (ML) for calciumrich rocks (including dolomites and anorthosites), mafic rocks, and gabbros for 5.25 < V P < 7.25. Figure 5 shows the comparison of the values so obtained with the measured data from the various authors. We considered the V P and V S computed by Faccenda et al. (2007), for four synthetic profiles from 0 to 22 km, taking into account the correction for temperature and pressure (F07 in Figure 5). We also considered the values of V P and V S from  (B12 in Figure 5) with a spread function <2.5, reducing the maximum depth considered to 12 km. Figure 5 shows that the seismic velocities (F07 and B12) fit well with the ones predicted by the ML. The results of the local earthquake tomography of Anselmi et al. (2011) in the Southern Alps confirm this observation. Viganò et al. (2013) results fit the ML only for depths shallower than 10 km in the eastern and south-western corners of the tomographic grid (around the Garda lake and in the eastern Southern-Alps). For greater depths and in other parts of the grid (in the Eastern Alps, and the area around GL), they better agree with the BRF. Behm (2009) computed a 3D model of the V S of the Eastern Alps and the Bohemian Massif from the tomographic inversion of stacked active-source data used to compute the V P model (Behm et al., 2007). Based on the information about V P and VS, Behm (2009) found that the Poisson's ratio in the Southern Alps and Dinarides between depths of 1 and 5 km is close to the values predicted by ML, while between 11 and 15 km, it is nearer to the values predicted by BRF. The Poisson's ratio in the Tauern Window (Figure 1) is close to the BRF values between 1 and 5 km depth, while it is lower than them between 6 and 15 km. In the Eastern Alps (excluding the Tauern Window) the Poisson's ratio values are grouped around the BRF line for the depth interval of 6-15 km, while at shallower depths the values are comprised between the BRF and ML lines. Therefore, we used the ML empirical relation in the eastern Southern Alps and the Dinarides from the topography up to 10 km of depth and the BRF law otherwise (inset in Figure 5).
For the bulk rock density, ρ, the information is also unevenly distributed. The laboratory measurements of Faccenda et al. (2007) add on the results of the gravimetric and seismogravimetric inversions in the area (Dal Moro et al., 1998;Ebbing et al., 2001Ebbing et al., , 2006Šumanovac et al., 2009;. After converting the V P into densities using the relation of Ludwig et al. (1970) as reported by Brocher (2005), we compared the resulting distribution with the V P and density values from the various measurements (Figure 6). We can see that Ludwig et al.'s (1970) Ludwig et al. (1970). The use of empirical relations to convert V P into V S and ρ introduced an additional source of uncertainties in our model. Brocher (2005) reported the average misfit of his dataset produced by different empirical relations, but he did not quantify the dispersion around the curves. We compared ML and BRF to the values in our dataset (Anselmi et al., 2011;Viganò et al., 2013) and we obtained an average model misfit of less or equal than 0.05 km/s, with a standard deviation of around 0.1 km/s for both relations.
We applied the conversion relation before the interpolation, and we analyzed the effect of this procedure. We interpolated V P from Anselmi et al. interpolated V P (V S 2). We compared V S 1 and V S 2 with V S 0: average misfits are comparable (-0.02 and 0.02 km/s), but the standard deviation of misfit of V S 2 (0.12 km/s) is twice of that of V S 1 (0.05 km/s), confirming greater effectiveness when empirical relationships are applied before the interpolation.
For the sediments, in the absence of laboratory measurements, we considered the properties reported by Giustiniani et al. (2008) for the Friulian plain. The values are in agreement with the measurements of the elastic moduli obtained from Vuan et al. (2011; the yellow area in Figure 2) and Tondi et al. (2019) for the Venetian and Po plains. Hence, we assigned these values to the grid points within the sedimentary layer ( Table 2). Figure 7 shows the 3D NAC multiparameter model: in this case, the V P is shown. Two vertical sections, one, in the X direction at y = 200 km (∼46.3 • N), the other in the Y direction at x = 240 km, i.e., (∼13.15 • E) enable to appreciate the vertical V P variations. The figure shows the Moho and an isosurface corresponding to V P = 6.4 km/s. Even if in stable continental regions commonly the lowermost part of the crust has V P greater or equal to 6.8 km/s (Mooney, 2010), in this area, the lower crust has a lower value of V P , being about 6.4 km/s at the transition between AD and EU (Bleibinhaus and Gebrande, 2006). Therefore, we interpreted the isosurface corresponding to a V P = 6.4 km/s as representative of the top of the lower crust, although NAC has not any vertical subdivision below the sedimentary layer.

MAIN FEATURES OF THE MODEL
We compared the input data with the final depth of the interfaces and the V P values by interpolating the values from NAC in the position of the data points (Figure 8). We found that there is no systematic over-or under-estimation of the data and that the distribution of the differences is compatible with the input data uncertainties. NAC shows a high variability of the geophysical parameters (  depths <10 km where V P ranges from 5.0 km/s to 6.5 km/s (with total uncertainties between 0.2 and 0.4 km/s) and estimated E has values from 50 to 100 10 9 Nm −2 (with total uncertainties between 5 and 20 10 9 Nm −2 ). The relative uncertainty is around 5% for V P , V S, and ρ, and between 5 and 20% for µ and E. Figure 9 shows more in detail the depth of the interfaces of NAC. The MO1 depth (9C) ranges from a minimum of 25 km in the southern part to a maximum of 55 km in the Alps. The maximum depth is found in the Central Eastern Alps, while two minima are present: one between the Garda Lake and the Adriatic Sea, and the other in north-eastern Slovenia and Austria. The western part of NAC, where the Adria Moho reaches the shallowest depths near the boundary with Eurasia, is characterized by a steep gradient of the Moho depth. In the eastern part, the transition is more gradual, as shown by the Alp01 profile  Figure 2). Figure 9 also shows the uncertainties of the final maps (σ TOT ). In general, the accuracy of the maps is variable, and it reflects the different distribution of the input data. The uncertainty due to errors in the input data (σ UNC ) shows low geographic variability and is comprised between 1 and 2 km.
The difference between the Moho depth of AD and EU at their boundary has its maximum value (up to 15 km) in the area close to the TRANSALP profile while it decreases from west to east, where it is comparable with the uncertainty of MO2. The differences between MO1 and MO2 are relevant only near the boundaries between the three Moho fragments, in particular in the western part of NAC between AD (MO2 shallower than MO1 up to 7 km) and EU (MO2 deeper than MO1 up to 10 km), and at the east of the triple junction, between PA and EU (differences up to 6 km) (Supplementary Figure S2).
In Figures 10A-C, two N-S profiles of the model (A-A' and B-B'), and one NE-SW oriented (C-C') enable the comparison of MO1 (black line), and MO2 (gray and blue lines, for EU and AD, and green for PA, respectively). Figure 10 enables to appreciate also the V P depth variations along the same profiles. The three sections show high variability both in the upper (0-20 km) and lower (z > 20 km) crust. The uncertainty (σ TOT ) is shown in Figures 10D-F: it is minimal for the central part of the model, with values <0.25 km/s. NAC's southern part is less constrained, especially for the section B-B' , with uncertainty even above 0.35 km/s. As inherent in the modeling procedure, the high velocity corresponds to an increase in ρ and the elastic parameters of the multiparametric model: in particular, zones with high velocity, ρ, and rigidity are present in correspondence of the Southern Alps and the External Dinarides front. An example is shown in Figure 11, presenting the vertical sections of the profile BB' of Figure 10 for V S (A), rho (B), µ (C), and E (D), below an exaggerated topographic profile. V S ranges between 3.0 and 4.0 km/s, with smaller values close to the front of the Southern Alps foothills. It is noteworthy, the thicker layer with V S around 3.8 km/s of Eurasia if compared with the southern part of the section. ρ shows a positive vertical gradient, with values ranging between 2.6 and 2.9 kg/m 3 . The elastic parameters  have a pattern similar to the one of V S , showing the maximum heterogeneity and the minimum values in correspondence of the Southern Alps' front. The value of the elastic parameters is greater in the southern part of the profile (AD) than in the northern part (EU). Figure 10 shows a similar pattern for V P , also comparing AD and PA (Figure 10C), even considering the error due to the uncertainty.

DISCUSSION
The building of NAC was driven by the need of a detailed model of the crust and its properties at the transition from the Alps to the Dinarides, where a vivid debate is ongoing on the relative plate movements, and, hence, the structure of the deep crust and of the upper mantle (Lippitsch et al., 2003;Kissling et al., 2006;  Mitterbauer et al., 2011;Handy et al., 2015;Kästle et al., 2020).
The model aims at serving as a basis for modeling the stress and strain fields in this critical area and interpreting the intense seismic activity of the region. To get new insight into the deep crustal structure, and on the processes having formed and deformed it, in the last decade several crustal models were developed at the continental scale of Europe, and in the last few years, seismic and seismological experiments are bringing new information. However, a detailed 3D model of the physical properties of the crust of the northern tip of AD is still lacking. The three most popular crustal models at European scale (EuCrust-07, EPCrust, and EUNAseis), with their uniform parameterization allowing a simple conversion in a numerical mesh, are an essential reference for constructing a regional model. Their scale and resolution yet make them unsuitable for studies in small areas, like the one we consider. EuCrust-07 has a resolution of 15' × 15' (or 0.25 • × 0.25 • ), EPCrust and EUNAseis of 0.5 • × 0.5 • , while a higher resolution characterizes the ESC Moho: 0.1 • × 0.1 • . Nonetheless, the authors of the latter, to avoid aliasing effects, transformed the original data points to 10 × 10 km block averages, and low-pass filtered them, using a cut-off length of 100 km, and a passing wavelength >200 km, to highlight the features at a continental scale. More detailed models, developed for seismological studies, cover only marginally our region (Vuan et al., 2011;Molinari et al., 2012Molinari et al., , 2015a. NAC fills this gap: it is a 3D model of the crust from the topographic surface to the Moho discontinuity. The variations in physical properties (seismic velocities, density ρ, Young modulus E, and shear modulus µ) are sampled on a grid with cells of 5 × 5 km horizontally, and 1 km vertically (Figure 7). A highly resolved Moho structure is essential for understanding the dynamic processes at plate boundaries (e.g., Carbonell et al., 2013, and references therein). With greater detail, however, the consistency of the information from different data sets, as, e.g., from different seismic experiments (refraction/wide-angle surveys, or receiver functions) becomes crucial. In most cases, the Moho depth from these two sets of data is consistent, although not always coincident, since the Moho depths from normal incidence are generally obtained from refraction derived models (Carbonell et al., 2013). The receiver function information is increasing with the time and the multiplication of experiments, and several works adopted different weights to combine it with the controlled seismic source data in the Moho definition (e.g., Di Stefano et al., 2011).

Selection of Input Dataset and Model Construction
NAC is mainly based on data from active seismic experiments that provided good coverage of most of the studied region (Figures 2, 3). We verified the consistency of their interpretation at the intersections (Supplementary Figure S1). The differences are mostly within the errors related to the original data acquisition and processing. The profiles with similar acquisition and processing characteristics (e.g., the profiles from ALP 2002 experiment), or interpreted by the same authors (e.g., the profiles 1 http://www.crs.inogs.it/bollettino/RSFVG from Scarascia and Cassinis, 1997) have a smaller difference between each other compared with the pairs that do not share these features (e.g., the TRANSALP and Alp02 profiles). The only way to remove such differences would be to reprocess all the profiles consistently. As observed by Brückl et al. (2007) based on the TRANSALP and Alp02 profiles, some differences in the velocity-depth functions could also be related to a degree of anisotropy in the velocity, and, therefore, are intrinsic in the different profile orientation. We integrated the dataset about the Moho depth with receiver function studies, overall for the parts of the study region not covered by the seismic profiles, like the south-western part of NAC (Figure 2).
We also considered the recent work of Hetényi et al. (2018), based on the analysis of the receiver functions along an NS transect at 13.33 • E, in the framework of the EASI complementary experiment of AlpArray. They applied different methods to obtain the Moho depth from the receiver functions inversion, obtaining a good agreement to the north of the Alps but highly scattered data in correspondence of the Eastern and Southern Alps, for the greater complexity of the crustal structure and possible reverberations. Therefore, Hetényi et al. (2018) chose to focus on the migrated receiver function profiles showing a very deep Moho (up to 75 km) on the Adriatic side. The proposed depth in this area is significantly different from the seismic profiles (Alp01, Alp02, and the three profiles from Scarascia and Cassinis, 1997), and other receiver function results in the same area (Bianchi et al., 2015) and also the results obtained by Spooner et al. (2019) with gravimetric methods. Hetényi et al. (2018) analyzed in detail the data below the Eastern Alps and reported that they do not find a sharp crust-mantle boundary but a broad zone with a velocity gradient, i.e., a clear Moho. However, the most plausible velocity contrast they found using the "west profile" (35-39 km) would be in agreement with the previous data, showing an Adria Moho shallower than the European one, and a Moho maximum depth <50 km. As recognized by the authors, the Moho depth range that can be obtained by their analysis is broad (between 50 and 70 km depth), and the model reaches mantle velocities over a broad depth range (48-65 km). Because of the non uniqueness of their results and the disagreement with other models in the same area, we decided to exclude from the dataset the results of Hetényi et al. (2018), although we considered the data in some of the following tests.
We choose to weight the information from different sources and to exclude part of the data based on their uncertainties. Therefore, uncertainty estimations are critical pieces of information for the construction of NAC. Unfortunately, only in a few cases (e.g., Behm et al., 2007), a careful analysis of uncertainty was done and reported for seismic experiments so that we had to rely on considerations on the applied techniques (Grad et al., 2009b). The estimations of uncertainty in the receiver function studies considered were performed in different ways, but they are based mainly on the bootstrap technique. Lombardi et al. (2008) also considered two different contributions to the total uncertainty, namely due to uncertainty of V P and finite bandwidth. They also concluded that the natural frequency content of the data is the major contributor to the total uncertainty of their results. Piana Agostinetti and Amato (2009) provided a qualitative classification of the reliability of the results.
Hence, the different approaches used in the literature to define the uncertainty make these estimates somehow subjective. To understand how the selection of data input and the weighting of uncertainty influenced our results, we performed three tests on MO1: (a) we interpolated all the data-points (also considering the data from Hetényi et al. (2018) without weights (MO1a), (b) we removed from the datasets the data from Hetényi et al. (2018) (MO1b), (c) we considered the weights based on the uncertainty (MO1c). The differences between the results of these tests and MO1, Z = MO1x − MO1, are reported in Supplementary  Figure S3. We made these tests on MO1 because the definition of the boundary between Eurasia and Adria is not compatible with the results of data from Hetényi et al. (2018). The main differences between MO1a and MO1 are along the profile of Hetényi et al. (2018) with Z < −20 km. Removing data from Hetényi et al. (2018) from the dataset, the differences of MO1b relative to MO1 are in the range of σ TOT , with the single exception of the SW corner of the NAC model, where the data coverage is poor. The inclusion of the weights due to data uncertainty further reduced the difference of MO1c relatively to MO1. As we can expect, the most unstable area coincides with the southwestern corner of NAC, where the results changed appreciably in the different tests, while the eastern part of NAC is stable also changing the selection of the receiver function data.
The velocity models of the considered seismic profiles have multiple layers for the crust (with the layer number different from one model to another) and a vertical velocity gradient within the layers. Some recent crustal models (e.g., EuCRUST-07; EPcrust) used a fixed number of layers for the whole area (sediments, upper crust, and lower crust) with horizontal variations of layer thickness and the velocity. This approach is useful and justified at the continental scale, but it is too simplified to show the real complexity of the studied region. Another problem with this approach is in the arbitrariness in associating the original layers of the seismic profiles to the upper crust or the lower crust since the most evident velocity contrasts are not always associated with the same interfaces. As we already said, the division adopted in stable continental regions, where most often the lowermost part of the crust has a V P greater or equal to 6.8 km/s, is not appropriate for our area. Therefore, we chose to adopt a single layer for the crust, in which, however, the seismic velocities, ρ, and the elastic parameters are variable, both laterally and vertically. In this way, NAC includes high-velocity bodies in the upper crust and low-velocity bodies in the lower crust.
An essential step in the construction of the NAC model was the estimation of model uncertainty on the depth of the interfaces and the geophysical parameters, taking account of the uncertainty of the input data (σ UNC ) and the effect of the interpolation process (σ INT ). The total uncertainty (σ TOT ) of the Moho depth (both MO1 and MO2) ranges from 2.5 to 6.5 km, which means that the relative uncertainty lies between 5 and 20% (Figure 9). The variability of the uncertainty is due mainly to the distribution of the input data. The map of the Moho uncertainty (Figures 9D,F) shows values mostly >5 km where the d min values are higher than 40 km, while the Moho uncertainty is comprised between 3.0 and 5.5 for smaller distances (see the map of d min in Figure 3).
In the construction of NAC, we had to deal with some problems typical of the spatial-data interpolation, like the inputdata selection, uncertainty and resolution evaluation, and the choice of the interpolation grid. Some of these issues could be better addressed by using a Bayesian approach where the solution is not a single best model but instead a probability distribution representing the full level of knowledge about model parameters (Bodin et al., 2012). The use of the transdimensional inversion (Sambridge et al., 2013) allows to a data-driven choice of parametrization, like the interpolation grid. A future improvement of NAC could be based on this approach.

Comparison Between MO1 and MO2
Another peculiarity of NAC is that in its construction we performed the estimate for the Moho discontinuity, both taking into account the hypothesis of a unique (MO1) as well of a fragmented Moho (MO2), to describe at best the Moho heterogeneity.
The knowledge on the relative position of the Adria and Eurasia Moho in the western part of NAC is mainly due to the results of the TRANSALP experiment, whereas eastwards, the interpretation of the seismic profiles is more questionable. The active seismic data at disposal from Scarascia and Cassinis (1997) and the ALP 2002 experiment  shows that the Adria Moho is shallower than the European one. However, as we said before, Hetényi et al. (2018) found a very deep Moho (up to 75 km) on the Adriatic side. It is worth noting that their boundary between AD and EU is further north than ours. If confirmed by other data, this would imply two oppositely dipping subduction zones (Király et al., 2018) and should lead to a re-analysis of the other datasets. The fragmentation of the Moho, leading to our MO2 model, was suggested based on the analysis of the seismic data of the CELEBRATION and ALP 2002 experiments (Behm et al., 2007) and supported by the elastic plate modeling of Brückl et al. (2010), gravimetric, geodetic, and seismicity pieces of evidence (Brückl, 2011). However, the existence of a distinct Pannonian fragment is not recognized by all the authors (e.g., Spada et al., 2013), even if the sharp boundaries between the strain velocities observed are in agreement with the triple junction between EU, AD, and PA (Völksen et al., 2018). The differences between our two Moho interfaces are small (Supplementary Figure S2) and they are more pronounced in the model areas well constrained by the data and negligible in the more controversial domains (see Figure 3). To the west, along the EU-AD plate boundary, MO2 is up to 7 km shallower than MO1 for AD, and deeper than MO1 up to 10 km for EU, whereas the differences at the east of the triple junction, between PA and EU reach a maximum of 6 km (Supplementary Figure S2). One could speculate that the choice of one or the other model is irrelevant for most of the applications. However, the analysis of the semivariograms can have a role in opting for one or the other model. The semivariograms for the different fragments are different in terms of the sill, range, and nugget; therefore, the stationarity assumption appears to be violated for MO1. AD shows a higher variability than the other fragments: the reason can be the greater heterogeneity of the used dataset or the greater complexity of Adria Moho with respect to EU and PA for the considered region. Moreover, the uncertainty σ TOT has a different distribution for MO1 ( Figure 9D) and MO2 ( Figure 9F): it is generally lower for MO2, except that in the model eastern and southern parts, and at the fragments boundaries. These differences depend on the semivariogram models used for MO1 and MO2 and on the distribution of input datasets. The lower values of MO2 are the effect of the lowest value of the nugget for the three MO2 fragments compared with the nugget of MO1. The higher values in the region less covered by the data, on the other side, depend on the lower ranges of AD, EU, and PA relative to the one of MO1. Another reason is that at some points along the boundary of the fragments, MO2 values are produced by an extrapolation, rather than by an interpolation inside the convex hull of the input datasets.

Comparison With Previous Models of the Moho
NAC shows a high variability in the Moho depth, especially in the AD domain, as it may be appreciated in Figures 9C,E, with a depth ranging from a minimum of 25 km between the Garda lake and the Adriatic sea to a maximum of 55 km in the Central Eastern Alps. Compared with other models at a broader scale (EuCrust-07; ESC Moho; EPCrust; Moho map of Spada et al., 2013), NAC is more detailed, with more pronounced minima and maxima. Although NAC does not contradict the principle of simplicity (Kissling, 1993), it evidences the differences in depth that can be interpreted in terms of different plates or fragments. To enable a correct comparison with the other models, it is necessary to use the same grid to resample a model to the grid of the other. When the grid spacing is appreciably different (EuCrust-07; EPCrust) we decided to resample NAC at the same spacing of the model with which we want to compare it, and vice versa: in this way it is possible to appreciate both the longwavelength variations, as well as the details of the model with the finer grid. When the grid spacing is similar (ESC Moho; Moho map of Spada et al., 2013), we used the grid of NAC. Depending on the model to be compared with, we considered MO1 or MO2, on a turn.
The long-wavelength features of the MO1 map of NAC are similar to the ESC Moho ( Figure 12A). The position of the maximum depth of MO1 and ESC Moho is well constrained by the profiles of Scarascia and Cassinis (1997) and TRANSALP. Superimposed to the long-wavelength, NAC shows short-wavelength variations, as the undulations in the high-Adriatic Sea. The other most striking difference occurs in the Southern Alps, where the differences between the Eurasia and Adria are more pronounced (about 10 km greater than in ESC Moho). This is reasonable because the inputs are similar, but ESC Moho was filtered passing wavelengths larger than 200 km. The EuCRUST-07 Moho (Figure 12B) locates the maximum depth of the Moho slightly to the north of the corresponding troughs of NAC and ESC Moho. The differences are relevant also in the eastern part of NAC, which benefits of the results of the more recent ALP 2002 seismic experiment. NAC and EuCRUST-07 also disagree on the depth in the south-western corner, while being very similar in correspondence of the Northern Adriatic Sea. To construct the EPCrust model in our area, Molinari and Morelli (2011) averaged information from EuCRUST-07, ESC Moho, and other regional models (Stehly et al., 2009;Molinari et al., 2010). The main difference with NAC ( Figure 12C) is, as for EuCRUST-07, the position of the maximum depth of the Moho, even if the shift is smaller than in Figure 12A. Finally, we compared MO1 with the Moho depth obtained by Spooner et al. (2019) (Figure 12D). The Moho of Spooner et al. (2019) is slightly shallower than MO1, but the two models are consistent: the only notable differences are in the south-western corner of NAC and in the central part of NAC, where MO1 shows a sharp transition between low and high depth values. The Moho of Spooner et al. (2019), aimed to fit the observed gravity in the Alpine region, was obtained trying to smooth the large vertical steps between defined plate domains, that are enhanced in models like MO2.
We compared our segmented Moho (MO2) with the Moho map of Spada et al. (2013) that also introduces the separation between EU and AD ( Figure 12E). The differences between the depths of the EU Moho are negligible, while the two models differ in the position of the boundary between the two plates. Moreover, Spada et al. (2013) did not take into account PA as a separated Moho surface. NAC's MO2 depth is systematically shallower than the one of Spada et al. (2013) in the AD domain (Central Southern Alps, eastern Po Plain, and Adriatic Sea), because of the constraints given by the data of Scarascia and Cassinis (1997) and of CROP experiment, not included in the other modeling.
Summarizing, except for the area less constrained by direct observations, the differences between our model and the other ones are sustained by the input dataset and are mainly related to the high resolution of NAC.
NAC allows choosing between the two models, one with a single Moho surface (MO1) and one with a fragmented Moho (MO2). The only limitation to the adoption of a fragmented Moho could be the more significant uncertainty of MO2 than MO1 in the eastern and southern part of the model (in AD) (Figure 9). However, if the chosen fragmentation is correct, the smaller uncertainty of MO1 is only apparent because we defined the Moho depth in a fragment with data from another plate. We are, therefore, in favor of MO2, although we provide also MO1, a smoother model that can be suited for, e.g., a gravimetric inversion.

Physical Properties of the Crust: Construction of the Model and Comparison With Independent Results
NAC is a multiparametric 3D model of V P , V S , ρ, µ, and E consistent over the whole area across North-Eastern Italy and Western Slovenia. Figures 10D-F allow to appreciate how the uncertainty of V P depends on the geographic distribution of the input data: it is minimal for the central part of the model, crossed by many seismic profiles, being higher in the southern part of the model, less covered by data (see Figure 3). We verified that the interpolation procedure that we used to obtain a uniform parameterization and resolution produced an unbiased model, compatible with the uncertainties of the input data.
A uniform resolution on the whole area means that the final model cannot reproduce all the fine details of the input data, as the high V P values (up to 6.8 km/s) reported by  in the upper crust of the Southern Alps. Nevertheless, compared with the model of Rossi et al. (2005), it is a much more representative model of the crustal heterogeneities in the area. We computed V S , ρ, µ, and E from V P values through empirical relationships. We preferred this solution to the alternative approach of constructing independent databases of V S and ρ and interpolating them. The reason is that the distribution and the resolution of the data are different for each parameter, and, therefore, the combination of their values needed to obtain elastic moduli can produce artifacts. We propagated the uncertainty of V P to the derived quantities, taking account also of the dispersion of V S and ρ values around the empirical relations. The relative uncertainty of V S and ρ is around 5%, while it is between 5 and 20% for µ and E.
As pointed out by Diaferia et al. (2019), applying an empirical relationship can be limiting, when the crust is heterogeneous, as shown by the high point dispersion (Figure 5). We tried to overcome this problem with the "regionalization" of the data, which resulted in the definition of two main domains, characterized by lithological differences. The results confirm the better fit of the ML relationship for the Southern Alps and the outer Dinarides from topography up to 10 km of depth, in agreement with the results of Behm (2009) about the relation between V S and V P in the different regions and the dominance of limestones and dolomitic limestones in the area . In the remaining of the model, where igneous and metamorphic rocks are dominant, and at greater depths, we used the BRF relationship. ML would be more appropriate for this region, if also here there is a pervasive fluidfilled cracks distribution at depth, according to the results of Diaferia et al. (2019) in the Apennines. The comparison with the available observation of V S values confirmed that we get greater effectiveness when empirical relationships to convert V P into V S , ρ, E, and µ are applied before the interpolation.
A possible drawback of this approach is that we assumed that the spatial variability of V P in the regions, where we applied the different empirical relations, provides valuable indications for the distribution of the most significant heterogeneities. Therefore, to check the validity of the procedure, we compared NAC with independent models. We compared the resulting V S of NAC with the values obtained independently by Kästle et al. (2018) through the joint inversion of ambient noise and earthquake phase velocity measurements in the whole Alpine area and by Guidarelli et al. (2017) through ambient noise tomography in an area including the central and the eastern parts of NAC ( Figure 13A). We re-interpolated the V S values of NAC in the same locations, considering only the crustal part of the other models. We found that the mean difference between Kästle et al. (2018) and NAC is -0.1 km/s while between Guidarelli et al. (2017) and NAC is about -0.2 km/s. In both cases, hence, the V S of NAC is slightly faster than the other ones but, considering that σ TOT of V S ranges from 0.1 to 0.25 km/s and the different methods applied, the comparison is satisfactory, especially for Kästle et al. (2018) values ( Figure 13B).
We also compared the density values of NAC with the ones from the crustal part of the model of Spooner et al. (2019). The values of NAC are systematically lower than the ones of Spooner et al. (2019) for depths >20 km ( Figure 13C). These differences and the high dispersion of density values around the Ludwig et al. (1970) relationship (Figure 6) suggest that the use of different empirical relations for different depths or a relation that takes account of the depth influence on the V P -density relation (as Christensen and Mooney, 1995) could give better results for the deepest part of our model. On the contrary, for depths ≤ 20 km, the differences between NAC and the model of Spooner et al. (2019) are minimal (mean of differences lower than 0.05 10 3 kg/m 3 ), confirming the validity of NAC in representing the physical properties in the region in this depth interval.

Relation Between Physical Properties of the Crust and Seismicity
NAC was built with the aim of modeling the stresses within the crust in relation to seismicity. The most seismically active region is the external front of the Southern Alps, and of the Dinarides. The seismic activity is concentrated at shallow depth (Viganò et al., 2015;Bressan et al., 2018a,b;Reiter et al., 2018). Also the Friuli, 1976 earthquake, the greatest event recorded in the instrumental period in this area, has a hypocentral depth of about 5-10 km (Slejko, 2018). , through the analysis of the spatial organization of seismicity in north-eastern Italy and western Slovenia, related the occurrence of seismicity in the shallowest crust to the mechanical rock heterogeneities. The interface separating media with different elastic moduli would be energetically favorable for the localization of fractures (Chatterjee and Mukhopadhyay, 2002;Bressan et al., 2018b). Spooner et al. (2019) compared the geographic distribution of large earthquakes in the Alps (M > 6) and found that their locations coincide with the upper-crustal density domain boundaries of their model, in agreement with the model of Ranalli (1992). To inquire on the relation between the seismicity pattern and the physical properties of the crust, we considered a vertical section of NAC crossing the region of the Southern Alps with the highest seismic activity and geodetic strain rate (Serpelloni et al., 2016), and well constrained by  (BB' in Figures 10B, 11). The pattern displayed is particularly complex at the front of the Southern Alps, where most of the seismic events are recorded, with a sharp transition from low to high values of seismic velocities, density, rigidity, and Young's modulus, in agreement with observations in other collisional systems (Ibarra et al., 2019). The analysis of the vertical sections AA' , BB' , and CC' supports the model of a strong and competent AD compared to softer EU and PA (e.g., Brückl et al., 2007Brückl et al., , 2010Marotta and Splendore, 2014).

CONCLUSION
We presented a new 3D model of the geophysical properties of the crust in northern Adria (NAC), by critically incorporating all the information available in the literature and extrapolating to the remaining parts of the region. NAC provides information on the depth geometry of the sub-sedimentary basement and the crust-mantle-boundary (Moho), and on the spatial variations both laterally and vertically of seismic velocities (V P , V S ), density (ρ), Young modulus (E), shear modulus (µ) within the subsedimentary crust. We quantified the uncertainties of interface depths and the physical properties by taking into account the data quality and the interpolation procedure.
NAC provides a higher resolution, especially in the Adriatic domain, than the existing crustal models at the scale of the European continent including this area. Differences in the Moho depth between NAC and the other models are sustained by the input dataset and are mainly related to the higher resolution of NAC.
We considered two different versions of the Moho surface: (i) one continuous surface (MO1) and (ii) three separate surfaces for the Adria microplate, Eurasia, and the Pannonian fragment (MO2). The differences between the two models are minimal, but, because the available data better sustain the solution of the fragmented crust MO2, we are in favor of MO2. We provide also MO1, a smoother model that can be suited for, e.g., a gravimetric inversion.
We computed V S , ρ, µ, and E from V P values by the application of empirical relations considering the regional distribution of different lithologies. We verified through a comparison with the available observation of V S values that the conversion of the original V P values should be done before interpolation. The comparison of NAC with independent models of V S and density confirms the validity of NAC in representing the physical properties in the region. At the same time, the results endorse the validity of the procedure followed, which can be applied in other contexts.
We developed NAC as a tool for a following modeling, aimed to analyze how the tectonic processes influence the intraplate stresses and the seismic activity. The model itself, however, already provides meaningful insights in this respect. The better fit obtained by the fragmented Moho, compared with the unique one, is an important reason for supporting this hypothesis. Similarly, the sharp transition from low to high values of seismic velocities, density, rigidity, and Young's modulus in the zone of the Southern Alps, where most of the seismic events are recorded, confirms the critical role of crustal strength gradients in conditioning deformation and seismicity patterns.
The model here presented can be precious for any other purpose requiring a detailed description of the physical properties consistent over the whole area, as studies on seismic wave propagation, geodynamic modeling, tomographic inversions. Moreover, it will be easy to improve NAC, including new information, when available. We also plan to test the possibility of adopting a data-driven grid to optimize the need for details with the reliability and availability of the pieces of information.

DATA AVAILABILITY STATEMENT
NAC digital model is available in the publicly accessible repository zenodo.org at link https://zenodo.org/record/ 3666708#.XnM51i2ZM6g.