Three-dimensional V_P and V_P/V_S models of the upper crust in the Friuli area (northeastern Italy)

Summary

3-D images of P velocity and P-to S-velocity ratio have been produced for the upper crust of the Friuli area (northeastern Italy) using local earthquake tomography. The data consist of 2565 P and 930 S arrival times of high quality. The best-fitting V_P and V_P/V_S 1-D models were computed before the 3-D inversion. V_P was measured on two rock samples representative of the investigated upper layers of the Friuli crust. The tomographic V_P model was used for modelling the gravity anomalies, by converting the velocity values into densities along three vertical cross-sections. The computed gravity anomalies were optimized with respect to the observed gravity anomalies. The crust investigated is characterized by sharp lateral and deep V_P and V_P/V_S anomalies that are associated with the complex geological structure. High V_P/V_S values are associated with highly fractured zones related to the main faulting pattern. The relocated seismicity is generally associated with sharp variations in the V_P/V_S anomalies. The V_P images show a high-velocity body below 6 km depth in the central part of the Friuli area, marked also by strong V_P/V_S heterogeneities, and this is interpreted as a tectonic wedge. Comparison with the distribution of earthquakes supports the hypothesis that the tectonic wedge controls most of the seismicity and can be considered to be the main seismogenic zone in the Friuli area.

earthquakes, gravity anomalies, seismotectonics, tomography

Introduction

3-D joint inversion for hypocentres and velocity structure (Thurber 1983; Eberhart-Phillips 1986; Evans et al. 1994) is a suitable method for investigating the regional tectonics and seismogenic characteristics of an area. The computed velocity images provide useful information on the relationship between superficial and deep tectonics, fault zones and the variation in seismic properties in the crust associated with the geological setting. Comparison with seismicity patterns provides the basis for a seismotectonic model explaining the relationships between the various structural units and the earthquake nucleation zones.

The upper crustal structure of the Friuli area (northeastern Italy) was first investigated by Amato et al. (1990) using local source tomography techniques. Their tomographic 3-D P-velocity model was characterized by strong lateral heterogeneities connected with a complex tectonic pattern. In particular, Amato et al. (1990) showed that most of the seismicity is localized within a southward-verging high-velocity body (V_P≥6.2 km s⁻¹), more than 5 km deep. Subsequently, Bressan et al. (1992), using the same tomographic inversion method but with a different set of data combined with gravimetric inversion, confirmed the previous results and interpreted the high-velocity body as the main seismogenetic zone of the Friuli area. All these investigations were based on P-wave analogue data.

In the present study we have computed 3-D V_P and V_P/V_S models, using better quality data (digitally acquired waveforms) from local earthquakes as recorded by the Friuli—Venezia Giulia local seismic network in order to investigate in detail the main tectonic features of the area.

Geological and Structural Setting

The geological units of the Friuli area (Fig.1) consist mainly of sedimentary rocks of Palaeozoic to Quaternary age. The Palaeozoic rocks are made up of limestones, terrigenous and volcanic deposits. Limestones and carbonate rocks characterize the Jurassic—Cretaceous period, while the prevalent Cenozoic and Quaternary deposits are flysch and molasse. The study area is characterized by a complex tectonic pattern (Fig.1), resulting from the superposition of several Cenozoic-age tectonic phases. The structural setting consists of two indented trapezoidal wedges, the outer containing the inner (Venturini 1991).

Figure 1

(a) Schematic geological map of the eastern Southern Alps. Symbols: (a) Hercynian low metamorphic basement (Ordovician—Carboniferous); (b) Palaeocarnic non-metamorphic chain (Upper Ordovician—Carboniferous) and Upper Carboniferous—Lower Permian covers; (c) Permo-Mesozoic, mainly carbonatic, successions; (d) flysch (Upper Maastrichtian—Middle Eocene) and molassic sequence (Miocene); (e) Quaternary alluvial deposits and moraines; (f) thrust; (g) subvertical fault. UD: Udine; TS: Trieste. Inset (b) Early syn-sedimentary faults reactivated during the Neoalpine compressions: PC—VB (Pieve di Cadore—Val Bordaglia fault); T—BC (Tramonti—But Chiarsò fault); D—I (Dogna—Idria fault); IL (Insubric Line); ML (Mojstrana—Ljubljana fault). The tectonic pattern outlined by these master faults consists of two indented tectonic wedges, named, respectively, inner and outer wedge. The grid of 3-D V_P and 3-D V_P/V_S inversions is shown (modified from Bressan et al. 1998).

The tectonic constraints of both wedges are NE—SW- and NW—SE-orientated palaeofault systems (PC—VB, T—BC, D—I, ML faults in Fig.1), active as strike-slip faults from Palaeozoic to Middle Eocene times, which were re-activated during the compressional Cenozoic phases. The Mesoalpine (Dinaric) NE—SW compression generated NW—SE-orientated thrusts during the Middle—Late Eocene. A Middle Miocene—earliest Pliocene N—S-orientated compression followed, causing severe shortening of the central area by means of south-verging thrusts and backthrusts. The palaeofault systems acted as strike-slip faults. Later, a Pliocene NW—SE-orientated compression generated NE—SW-trending thrusts and folds, re-activating the NE—SW-orientated palaeofaults as inverse faults. Neotectonic activity was intense mostly in the southernmost belts of the inner and outer wedges. A schematic geological cross-section of the central part of the study area is shown in Fig.2.

Figure 2

Interpreted geological cross-section 7 N—S of Fig.1. Symbols: (a) Quaternary and Cenozoic units; (b) Mesozoic units; (c) Palaeozoic units; (d) fault.

The present state of tectonic stress, as inferred from the focal mechanisms of earthquakes (Bressan et al. 1998), appears to be significantly conditioned by the geometry of the indented wedges. It is characterized by a general compressional state of stress with a strike-slip stress state localized in the eastern part. The maximum compressional axis of stress is orientated from NW—SE in the western part to N—S in the eastern part of the area.

The seismicity is strongly focused in the central part of the Friuli area, with the maximum density of earthquakes occurring between 7 and 11 km depth. This area has been affected in the past by destructive earthquakes (Slejko et al. 1989) such as the 1348 event with epicentral intensity I₀ = IX MCS (Mercalli—Cancani—Sieberg scale), the 1511 event (I₀ = IX MCS), the 1928 event (I₀ = IX MCS) and the 1976 event (I₀ = X MCS, M_L = 6.4).

3-D Tomography

Inversion method and data set

The 3-D V_P and V_P/V_S tomographic images of the upper crust are obtained using Thurber’s (1983) method of joint inversion for hypocentres and 3-D velocity structure from local earthquakes. The computer code simulps (Evans et al. 1994) is used. A 1-D starting velocity model is assigned to the nodes of a 3-D grid and the velocity value at any other point is obtained by linear interpolation between the nodes. The traveltime from hypocentre to station is calculated by a pseudobending method (Um & Thurber 1987). The solution is then obtained by an iterative procedure, solving for hypocentre location and calculating the velocity of the medium, with a damped least-squares approach.

The 3-D V_P and V_P/V_S inversions are based on the arrival times of P and S waves of local earthquakes, recorded by the Friuli—Venezia Giulia seismic network (Fig.3). The seismic network consists of 15 short-period seismometers, with a digital acquisition system and a trigger algorithm for event detection. Data are sent via a telemetry radio link to the acquisition centre at UD (Udine, Fig.1). All stations are equipped with a vertical 1 Hz seismometer, except stations CAE, DRE and ZOU which are equipped with three-component 1 Hz seismometers. Since June 1995, station BAD has been equipped with a 3-D seismometer. From 1988 to April 1994 the sample rate was 360 sps on vertical sensors and 120 sps on each of the 3-D sensors, with a dynamic range extending up to 96 dB and 10-bit digitizing. The network was upgraded in April 1994 to include digital signal processing with a 16-bit digitizer, 120 dB of dynamic range and a sampling rate of 62.5 sps.

Figure 3

Location of the 3-D inversion grid: the triangles indicate the seismic stations; black dots, the earthquakes used for the V_P inversion; and white dots, those used for the V_P/V_S inversion. South to north grid nodes are marked as letters A, B, C, D, E, F, G, H, corresponding to Y (km) distance. West to east grid nodes are marked as numbers 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, corresponding to X (km) distance. (See text).

The area investigated is represented digitally by a grid extending 90 km in the west—east direction (west to east grid nodes at X = −50, −35, −25, −15, −7, 0, 7, 15, 25, 40 km) and 55 km in the south—north direction (south to north grid nodes at Y = −35, −20, −10, −5, 0, 5, 10, 20 km). The depth grid spacing is Z = 0, 2, 4, 6, 8, 12, 15 km. The coordinates of the grid centre are latitude 46°20́N and longitude 13°05́E. The horizontal mesh of the grid is finer in the central part to give an approximately homogeneous ray coverage within the blocks. The data set consists of local earthquakes with local magnitude ranging from 1.3 to 4.3, recorded by the seismic network from 1988 to 1997, located with a HYPO71 program (Lee & Lahr 1975), with GAP (largest azimuthal separation between stations) ≤180°. The 3-D V_P inversion was obtained using 2565 P-wave arrivals from 224 events, with an estimated picking accuracy of ±0.05 s. The data set for the 3-D V_P/V_S model consists of 930 S-wave arrival times from 136 earthquakes, most of them recorded with the upgraded acquisition system, which allows a picking accuracy of ±0.1 s for 3-D seismometers and ±0.15 s for vertical seismometers. The earthquakes for the V_P/V_S model were carefully selected by visual inspection to ensure a significant number of S arrival times.

Inversion procedure

Given that the S-wave data are of poorer quality and fewer in number than the P-wave data, we started by constructing a reliable 3-D V_P model and then inverting for V_P/V_S values, based on the results of the 3-D V_P model.

The reliability of the results obtained with linear tomographic inversion depends on the initial reference model. Kissling et al. (1994) showed that an inappropriate initial reference model could significantly affect the quality of the tomographic images. We therefore followed the Kissling et al. (1994) approach by obtaining an initial reference model (the minimum 1-D model), using the velest program (Kissling et al. 1995).

An a prioriP-velocity model of the Friuli area with constant-velocity layers was compiled using existing refraction seismic data (Pagini 1995). The minimum 1-D V_P model was then computed with the VELEST program by inverting for hypocentres every iteration, and inverting station delays and velocity values every second iteration until the total RMS value (root mean square misfit of traveltime residuals) reduced significantly and stabilized. Table 1 shows the a priori and the calculated minimum 1-D models. The layer at negative 3 km depth has been included to account for the Earth’s topography. The damping value for the 3-D inversion is 10, selected following the empirical approach of Eberhart-Phillips (1986) to give both low data variance and low standard error. The effect of topography in the 3-D inversions has been considered by including the elevation of the stations.

Table 1

1-D V_Pa priori model and the minimum 1-D V_P model computed with the VELEST program. The first layer (−3.0 to 0.1 km depth) includes all the seismic stations.

The 3-D RMS residual is 0.203 s. Fig.4(a) shows the V_P images at six depth slices. The diagonal elements of the resolution matrix and the standard error are shown in Figs 4(b) and (c), respectively. Since the resolution is also controlled by the ray density, the best resolution is found in the central part of the grid at depths from 2 to 6 km, where the station coverage is denser and the number of earthquakes is larger. The poor resolution in other sectors of the grid is explained by the poor density of earthquake ray paths.

Figure 4

Results of the 3-D V_P inversion at six depth slices (0, 2, 4, 6, 8 and 10 km). The X (km) and Y (km) distances are marked (see Fig.3). (a) P-velocity images; (b) diagonal elements of the resolution matrix; (c) standard error.

The V_P/V_S minimum 1-D model; that is, the initial reference model for 3-D V_P/V_S inversion, was also computed with the velest program (Kissling et al. 1994). The a priori 1-D model consists of the average V_P computed for each layer from the calculated 3-D V_P model, and V_S values taken from Riggio et al. (1987) and Fah et al. (1993). The earthquakes were relocated by fixing the V_P values and inverting for the shear wave velocity only. The a priori V_S model and the calculated minimum V_P/V_S 1-D models are shown in Table 2. The damping value for the 3-D V_P/V_S inversion, based on the Eberhart-Phillips (1986) method, is 5. The 3-D V_P/V_S data were inverted by fixing the V_P values as obtained from the 3-D V_P model and allowing hypocentre relocation. The 3-D RMS residual is 0.274 s. The plan-view maps of the V_P/V_S anomalies are shown in Fig.5(a). Fig.5(b) and (c) are, respectively, the images of the diagonal-elements resolution and the standard error. The central part of the area shows higher resolution at depths from 2 to 6 km, as expected from the station coverage and distribution of earthquakes.

Table 2

1-D V_Sa priori model, the minimum 1-D V_S model and the 1-D V_P/V_S model computed with the VELEST program.

Figure 5

Results of the 3-D V_P/V_S inversion at six depth slices (0, 2, 4, 6, 8 and 10 km). The X (km) and Y (km) distances are marked (see Fig.3). (a) P- and S-velocity ratio images; (b) diagonal elements of the resolution matrix; (c) standard error.

Since most of the S arrivals are taken from vertical sensors, the V_P/V_S inversion is subject to larger errors than the V_P-only inversion. The picking of S waves from vertical sensors is less accurate than for the P arrivals because they are partially obscured by the P-wave signal. However, a more serious problem in identifying shear waves from vertical instruments (Thurber & Atre 1993) arises from S-to-P converted waves, which arrive before the true S wave and which can be wrongly identified as the direct signal. This can result in artificially low values of V_P/V_S, especially in areas with unconsolidated sediments overlying rock units with sharp velocity contrasts. This effect is negligible in the present study because only two vertical stations (TLI, UDI) are located on unconsolidated sediments (Fig.3), and the number of S arrival times from these sites used in the V_P/V_S inversion is 24, a small fraction of the total employed phases (930).

The calculated 3-D V_P/V_S model (Fig.5a) shows values decreasing with increasing depth, a reasonable trend consistent with the closing of microcracks under increasing confining pressure.

Resolution analysis

Following the approach of Zhao et al. (1992), we analysed how the true velocity images were reconstructed in the inverted velocity images using the checkerboard resolution test (CRT) and restoring resolution test (RRT).

The goal of the CRT is to investigate the adequacy of the ray coverage and the resolution. The checkerboard model was constructed by assigning to the 3-D grid nodes positive and negative fractional V_P (6 per cent) and V_P/V_S (3.5 per cent) perturbations. Figs 6(a) and (b) show V_P and V_P/V_S fractional perturbations in per cent, obtained with the 3-D inversion. Figs 7(a) and (b) are, respectively, the images of V_P and V_P/V_S synthetic checkerboard inversions. The V_P checkerboard pattern is well reconstructed in the central part of the area for layers at 2, 4, 6 and 8 km depth, with the best reconstructed patterns at 4 and 6 km depth. The resolution decreases at 10 km depth because of the low number of quakes, with consequent relatively poor sampling of nodes by ray paths. The checkerboard pattern is not correctly reconstructed at 0 km depth because of uneven sampling of ray paths. The resolution of the V_P/V_S images is lower than that of the V_P images, and the best reconstructed patterns are at 4 and 6 km depth.

Figure 6

Fractional (a) V_P and (b) V_P/V_S perturbations (in per cent), obtained with the tomographic inversion at six depth slices (0, 2, 4, 6 and 8 km). The X (km) and Y (km) distances are marked (see Fig.3).

Figure 7

Results of the checkerboard resolution test (CRT) for (a) V_P images and (b) V_P/V_S images. The X (km) and Y (km) distances are marked (see Fig.3).

The goal of the RRT is to check how errors in the data set influence the inverted images and the relocation of the earthquakes. The 3-D tomographic images obtained were considered as the synthetic model of the RRT. Random errors with a normal distribution, as in the real data set, were added to the calculated synthetic arrival times. A restored image was constructed by inverting the synthetic data. The RRT V_P images (Fig.8a) show that the tomographic images are well restored for layers at 2, 4, 6 and 8 km depth, while the 3-D inverted V_P/V_S images are less well restored (Fig.8b). The best restored V_P/V_S tomographic images from the RRT are in the central parts of the layers at 4 and 6 km depth. The average deviation from the initial locations of the hypocentre locations obtained with the RRT is 0.3 km.

Figure 8

Results of the restore resolution test (RRT) for (a) V_P images and (b) V_P/V_S images. The X (km) and Y (km) distances are marked (see Fig.3).

Discussion

Tomographic results

Based on the 3-D V_P/V_S model, we relocated 415 events that occurred within the study area during the period 1993–1996, in order to give an improved representation of the seismicity pattern for comparison with the calculated tomographic model.

The distribution of V_P and V_P/V_S anomalies is characterized by marked lateral and depth variations (Figs 9a to d), reflecting structural heterogeneities (Fig.1). Low P velocities (V_P = 5.4–5.8 km s⁻¹) and high V_P/V_S values (1.82–1.85) are related to superficial molasse and flysch deposits in the southern sector of the investigated area. Mesozoic limestones and dolomitic rocks, which constitute the main portion of the investigated crust, are characterized by a wider range of P-wave velocities (5.9–6.6 km s^‒1) and V_P/V_S values (1.78–1.88). V_P values between 6.0 and 6.5 km s^‒1 and V_P/V_S values in the range 1.75–1.82 below 5 km depth in the northern sector could be related to the Palaeozoic units (terrigenous sediments, limestone deposits, volcanic and low-grade metamorphic rocks). According to Cati et al. (1987), the pattern of magnetic anomalies indicates a crystalline basement at a depth of nearly 10 km, affected by south-verging thrusting and detached from the overlying Mesozoic cover. The resolution of the tomographic V_P and V_P/V_S models at depths greater than 10 km is poor and does not allow detailed images to be constructed. However, this zone is characterized by V_P values of around 6.3–6.5 km s^‒1, and V_P/V_S ratios range from 1.75 to 1.80.

Figure 9

Vertical cross-sections of the 3-D V_P and 3-D V_P/V_S images. (a) section 5; (b) section 6; (c) section 7; (d) section E. The V_P values are shown as contour lines. The V_P/V_S values are plotted in graded colours in (a) to (c), and in grey tones in (d). Diamonds indicate the positions of the relocated 415 earthquakes as calculated from the 3-D V_P/V_S model. The focal mechanisms of the strongest earthquakes of the 1976 sequence are shown in section 7: (1) May 6, M_L = 6.4; (2) September 11, M_L = 5.1; (3) September 15 M_L = 5.8; (4) September 15, M_L = 4.7; (5) September 15, M_L = 6.1. On section E, the traces of the main faults are shown: (a) Tramonti—But Chiarsò fault; (b) (c) Tagliamento—Osoppo fault system; (d) Dogna—Idria fault. The locations of the cross-sections are shown in Figs 1 and 3.

As expected, V_P/V_S values decrease with depth since cracks close with increasing pressure (O’Connell & Budiansky 1974). The near-surface high V_P/V_S values could result from the presence of water-saturated cracks, while the anomalously high V_P/V_S values (1.84–1.90) at depths greater than 3 km are attributed to highly fractured zones.

The most important feature revealed by the tomographic images is a high-velocity zone (V_P≥6.2 km s^‒1), bounded by sharp lateral and vertical velocity variations and located at about 6 km depth in the central part of the inner tectonic wedge (Figs 9a to d). Its base is not clearly resolved, but the boundaries of this zone are characterized by sharp V_P/V_S variations. The shape of the high-velocity zone is not regular. This is particularly evident in Fig.9(d), where the high-velocity body deepens towards the east. Most of the relocated earthquakes occur within or near the high-velocity body and along the high-V_P/V_S anomaly gradients (Figs 9a to d). We interpret this tectonic structure to be the result of the severe crustal shortening caused by the Middle Miocene—earliest Pliocene N—S-orientated compression (Venturini 1991). The V_P/V_S heterogeneities probably result from deformation along faults affecting the high-velocity body. The stronger shocks of the 1976 sequence (Barbano et al. 1985), the most important to occur within the area during the last century, are also shown in Fig.(c). All the earthquakes characterized by thrust focal mechanisms are located at the southern border of the high-velocity bulge.

The E—W-orientated tomographic images (Fig.9d) show evidence of the western border of the inner tectonic wedge, marked by the Tramonti—But Chiarsò fault. The Tagliamento—Osoppo strike-slip faults are marked by sharp lateral V_P and V_P/V_S changes in the central part of the inner tectonic wedge and clearly affect the geometry of the high-velocity zone. There is no clear evidence of correlation between the Idria—Dogna fault and the tomographic images because of poor ray coverage in this zone. The earthquake clustering that is noticeable in the western part of Fig.(d) consists largely of the swarm that occurred from 1996 January to June, with main shock local magnitude 4.3.

Excluding effects due to high temperatures, which are not expected from the geothermal environment of the area (average heat flow 50–60 mW m^‒2, Cataldi et al. 1995), the high V_P/V_S values at depth are attributable to a high degree of fracturing and high fluid pressure. Generally, we relate the sharp lateral V_P/V_S variations to the effects of faulting. The fluid overpressure can cause chemical effects, such as stress corrosion and pressure solution, with subsequent concentration and enhancement of local stresses, favouring the nucleation of earthquakes (Zhao & Negishi 1998).

Generally, the overall V_P and V_P/V_S pattern reflects the complex tectonic-structural setting, resulting from the superposition of three main tectonic phases (Venturini 1991). Each tectonic phase, characterized by different orientations of the principal axes of stress, inherited the deformations of the previous phase and re-activated the main regional faults. This has resulted in a complex 3-D deformation pattern, with discontinuous blocks and bulges marked by rapid spatial variations in V_P and V_P/V_S values. The high-velocity body is taken to represent the more brittle and stronger parts of the seismically active layers and is considered to be the main seismogenetic zone.

Experimental V_P measurements

For comparison with the tomographic results, laboratory measurements of P velocity were performed on samples collected in the ‘Dolomia Principale’ (dolomitic rock) and the ‘Calcare of Dachstein’ (limestone) units. These lithologies are representative of the Mesozoic units that constitute most of the crustal rocks in the study area.

The velocity of compressional elastic waves was measured at confining pressures of up to 300 MPa on cores of 26 mm in diameter and 35–55 mm in length, using the pulse transmission technique (Birch 1960). Sample preparation and details of the technique used are the same as in Burlini et al. (1998). Bulk density was determined from weights. The confining pressure/velocity relations are shown in Fig.10. The non-linear part of the curve is generally interpreted as being due to crack and pore closure (Birch 1961), while the linear part (at high pressure) reflects the intrinsic seismic properties of the rocks, that is the crack-free matrix properties, which correspond to the maximum velocities possible for the rock type under consideration. The maximum velocity at high confining pressure is 6.0 km s^‒1 for the limestone, and 6.9 km s^‒1 for the dolomitic rock.

Figure 10

Confining pressure/velocity relations of two cores of Friuli limestone cut parallel to layering (two mutually orthogonal directions, 1 and 2) and one core of Friuli dolomia cut normal to layering, as determined by laboratory measurement. Open symbols are increasing pressure measurements, full symbols are decreasing pressure measurements.

We computed the distribution of V_P, density and temperature with depth for a simplified lithological cross-section derived from Castellarin et al. (1979) (Fig.11). We calculated geotherms using a finite difference algorithm (1-D model) with the parameters listed in Table 3, assuming a temperature of 273 K and a heat flow of 50 mW m^‒2 at the surface. For the lower crust and the upper mantle we also used a second temperature derivative. The temperature obtained at each depth level was used to calculate the effective density in dry conditions, using the pressure and temperature derivatives, and thus the pressure at depth. These pressures and temperatures were used to calculate the seismic velocity at each depth level based on the experimentally determined velocities and their pressure derivatives. For the temperature correction of V_P we took derivatives as reported in the literature (Carmichael 1989).

Figure 11

V_P, density and temperature with depth calculated for the Friuli section. V_P is in km s^‒1, density in Mg m^‒3 and temperature in hundreds of °C. A simplified lithological cross-section derived from Castellarin et al. (1979) is considered: mainly limestones (from surface to 5 km depth); mainly dolomitic rocks (5–8 km depth); Palaeozoic rocks (8–11 km depth); quartzitic basement (11–12 km depth); upper crust (12–15 km depth); middle crust (15–25 km depth); lower crust (25–31 km depth); upper mantle.

Table 3

Thermal conductivity and its pressure and temperature derivatives, heat production, density and its pressure and temperature derivatives used to calculate the geotherm. These parameters were taken from a variety of sources (Haenel et al. 1988; Fountain et al. 1987 and references therein), and applied to the lithological section of Friuli with the same considerations as reported in Burlini (1992).

Since the 3-D V_P inversion is best resolved from the surface to 8 km depth in the central part of the study area, we examined the V_P/depth profile evaluated from laboratory data in this depth range for comparison with tomographic data. The V_P values (about 6 km s^‒1) from the surface to 5 km depth correspond to the average V_P values obtained with the 3-D inversion for this depth. The step at 5 km depth and the increase in V_P up to 6.8 km s^‒1 are in agreement with the presence of the high-velocity body, supporting the interpretation that the high-velocity body, or at least its uppermost layers, is mainly composed of dolomitic rocks. The calculated density of the limestones ranges from 2.64 to 2.69 Mg m^‒3, while the computed density of the dolomitic rocks varies from 2.79 to 2.82 Mg m^‒3.

Gravity modelling

The tomographic images were used for modelling Bouguer gravity anomalies, provided as digital data on a 3 × 3 km grid (Fig.12) by Cassano et al. (1989). 2-D gravity models were constructed along tomographic sections 5, 7 and E of Fig.9. The corresponding cross-section gravity models are given in Figs 13(a), (b) and (c), respectively. The P-velocity values were resampled along each section using a linear interpolation algorithm (De Boor 1978) with meshes spaced 2.5 km horizontally and 1 km in depth. The P velocities were converted into densities using the velocity—density relationship of Zelt & Smith (1992). The arithmetic mean density was calculated from the four density values of each mesh. The gravity response was calculated with the Talwani algorithm (Telford et al. 1990), using 2.85 Mg m^‒3 as background density. The model also includes layers representing the upper crust extended to the top of the lower crust (2.8 Mg m^‒3),the lower crust (2.91 Mg m^‒3) and the upper mantle (3.10 Mg m^‒3) as derived from Italian Explosion Seismology Group and the Institute of Geophysics, ETH Zurich. (1981).

Figure 12

Bouguer anomaly contour map (5 mgal contour interval, reduction density 2.4 Mg m^‒3). The grey box indicates the tomographic study area. The modelled gravity lines along sections 5, 7 and E (see Figs 1 and 9) are shown.

Figure 13

Gravity models of sections 5 (a), 7 (b) and E (c). Upper panel: anomaly profiles. Lower panel: density model with values in Mg m^‒3. The contour values of density variation of the optimized models are shown. In the gravity modelling a layer with density 2.8 Mg m^‒3 is included from 12?km depth to the lower crustal boundary. The lower crustal boundary is modelled with a north-dipping plane, located at 25 and 19 km depth from the surface, respectively, at the northern and southern edges of the tomographic grid. The density assigned to this layer is 2.91 Mg m^‒3. The upper mantle boundary is marked with a north-dipping plane located 40 and 35 km depth from the surface, respectively, at the northern and southern edges of the tomographic grid. The density is 3.10 Mg m^‒3.

The initial computed gravity anomalies and the observed N—S sections (Figs 13a and b) display a high degree of correlation. The mean residuals of the initial computed gravity anomalies are 4.2 mgal (section 5) and 3.8 mgal (section 7). Optimized density models, obtained using the least-squares inversion (More 1977), show mean residuals of 0.3 mgal for both sections. The maximum negative variations (‒0.08 Mg m^‒3 and ‒0.06 Mg m^‒3) of the optimized density model with respect to the initial density model are observed in the superficial layers of the distance range ‒25 to ‒15 km for both N—S sections. This disagreement could be due to the poor sampling of the shallow structure, which is not fully resolved by the tomography. The minimum density values (2.55 Mg m^‒3) of these layers are due to the presence of the shallow flysch and molasse sediments (maximum thickness 2.5 km) of the Friuli plain (Fig.2). The maximum positive variations (0.06 Mg m^‒3) are found in the shallow part of the modelled section 5, between 0 and 10 km distance, and are associated with density value of about 2.75 Mg m^‒3, appropriate to the Mesozoic rocks. Fig.13(c) shows the gravity response modelled along the E—W-orientated tomographic section E. The mean residual of the initial model is 4.5 mgal, and the optimized density model shows a mean residual of 0.8 mgal. The most noticeable difference between the computed and the observed gravity anomalies is found in the 30–40 km distance range, where the tomographic images show poor resolution because of poor sampling of the seismic ray paths. The other minor gravity discrepancies observed in the gravity profiles may well be due to local variations in the velocity—density relationship used. The prevalent geological unit consists of alternating Mesozoic limestones and dolomitic rocks (Fig.2). Considering the results from laboratory measurements and the distribution of V_P/V_S anomalies, our interpretation is that the strong density variations seen in the gravity models for the crustal layers investigated is due to fracturing and a change in lithology from limestones to dolomitic rocks.

Fig.14 shows a map of the horizontal gradient modulus computed from the Bouguer anomaly data. The high-gradient zone (>5 mgal km^‒1) corresponds to the boundary between the Friuli plain and the South-Alpine thrust belts. This zone is orientated NE—SW in the western part of the study area and E—W in the central part, marking the limits of the high-velocity body revealed by the tomographic inversion.

Figure 14

Modulus of the Bouguer anomaly horizontal gradient. The grey areas indicate gradients greater than 5 mgal km^‒1. The bold line indicates the interpreted limit of the high-velocity wedge. The modelled gravity lines along sections 5, 7 and E (see Figs 1 and 9) are shown.

Conclusions

The 3-D V_P and V_P/V_S images of the Friuli area show that the upper crust is characterized by marked heterogeneities that are related to the complex tectonic pattern. The V_P/V_S anomalies, and especially the sharp lateral V_P/V_S variations, appear to be related to the degree of fracturing, and hence to faulting geometries. The relocated seismicity occurs mainly along the high-V_P/V_S anomaly gradients. The V_P anomalies show that a high-velocity body is present in the central part of the area, at about 6 km depth, bounded by gradients in the V_P/V_S values. This zone is interpreted as a tectonic wedge, accommodating large amount of crustal shortening associated with Alpine tectonic phases. Most of the relocated seismicity and the stronger shocks of the 1976 sequence are located inside and around the tectonic wedge, which we consider to be the main seismogenic zone of the area. In accordance with laboratory V_P measurements performed on representative rock samples of the sedimentary crust, we relate the increase in V_P at the tectonic wedge to the presence of prevailing dolomitic rocks. The observed gravity anomalies show trends correlated with those computed using a density model derived from the 3-D V_P values. The final optimized density model is characterized by density variations which we ascribe to fracturing and lithological variations in the upper sedimentary layers.

Acknowledgements

We thank E. Kissling for kindly providing the velest program, D. Eberhart-Phillips for her suggestions, and A. Michelini and C. Venturini for helpful discussions. Thanks are also due to P. L. Bragato and M. Sedmach for elaboration of Figs 2, 9, 10 and 11.

The local seismic network is managed by the Dip. Centro Ricerche Sismologiche of the Osservatorio Geofisico Sperimentale with financial support from the Regione Friuli-Venezia Giulia.