High-resolution seismic reflection survey crossing the Insubric Line into the Ivrea-Verbano Zone: Novel approaches for interpreting the seismic response of steeply dipping structures

A high-resolution seismic reflection survey has been conducted across the Insubric Line from the Sesia Zone into the Ivrea-Verbano Zone (IVZ), where a remarkably complete cross-section of lower continental crust is exposed. The survey was carried out in preparation for the DIVE (Drilling the Ivrea-Verbano zonE) project, which was recently approved by the International Continental Scientific Drilling Program (ICDP). DIVE aims to gain new insights into the characteristics of the lower continental crust through targeted drilling, sampling, and borehole logging. A key borehole is planned near the Insubric Line at Balmuccia, where the deepest parts of the lower continental crust are exposed. As such, the primary objective of this seismic survey was to explore whether the sub-vertical structures prevailing at the surface can be expected to continue at depth or whether there are any indications for major flattening or fault-related offsets. Correspondingly, the acquisition and processing of the seismic reflection data were geared towards revealing weak backscattered events from local heterogeneities associated with the prevailing sub-vertical structural grain. The migrated sections, contain coherent back- scattered events to a depth of ~1 km, which form numerous short lineaments that seem to align sub-vertically. To substantiate this observation, we have generated synthetic seismic reflection surveys for canonical models of sub-vertical structures associated with Gaussian- and binary-distributed heterogeneities. Both the observed and synthetic seismic data were then subjected to energy-based attribute analysis as well as geostatistical estimations of the structural aspect ratios and the associated dips. The results of these quantitative interpretation approaches are indicative of the overall consistency between the synthetic and the observed seismic data and, hence, support the original qualitative interpretation of the latter in that the sub-vertical structural grain evident at the surface seems to prevail throughout the imaged part of the upper crust.


Introduction
The Insubric Line corresponds to the western end of the Periadriatic fault system delineating the boundary between the European and Adriatic plates, which can be followed over more than 700 km from Slovenia to the southwest of Torino (e.g., Schmid et al., 1989;Handy et al., 2015). Amongst the exposed paleo-plate boundaries, the Insubric Line is unique as it represents a particularly well-preserved suture zone documenting the late continental collision stage of the Alpine orogeny (e.g., Schmid et al., 1987). As such, the Insubric Line separates the Austroalpine gneiss units of the Sesia Zone to the west from Adriatic lower crustal rocks exposed in the Ivrea-Verbano Zone (IVZ) to the east.
The IVZ is associated with a very prominent positive gravity anomaly (e. g., Berckhemer, 1968;Kissling, 1984;Scarponi et al., 2020) and an associated seismic high-velocity anomaly (e.g., Diehl et al., 2009;Lu et al., 2018), whose origins extend to shallow crustal levels. This socalled Ivrea Geophysical Body (IGB) is assumed to represent a sliver of Adriatic lower crust and upper mantle, which was "upwarped" in the course of the continental collision process (e.g., Schmid et al., 2017). The IVZ thus corresponds to a partial exposure of the IGB.
The geology of the IVZ has been extensively studied, as it is widely regarded to represent arguably the most complete exhumed crosssection of the lower continental crust and its transition into the uppermost mantle (e.g., Fountain, 1976;Schmid et al., 1996;Quick et al., 2003;Brack et al., 2010) (Fig. 1). As a consequence, the IVZ is a prime locality to gain insights into the composition, structure, and evolution of this critically important, but still largely enigmatic part of the continental lithosphere. These questions will be addressed in a unified and focused manner in the framework of an international research project called DIVE (Drilling the Ivrea-Verbano zonE) (Pistone et al., 2017), which recently has been approved by the International Continental Scientific Drilling Program (ICDP). A ~4-km-deep borehole is planned just east of the Insubric Line near Balmuccia, where the deepest parts of the lower continental crust have been exhumed. To locate and orient this borehole, and to develop an effective drilling, sampling, and logging strategy, the targeted subsurface region needs to be geophysically characterized. This is generally best achieved through seismic reflection surveys, which are indeed widely regarded as a conditio sine qua non for drilling campaigns in general and ICDP projects in particular (e.g., Demirel-Schlueter et al., 2005;Juhlin et al., 2010;Simon and Buske, 2017). To this end, we have carried out high-resolution seismic reflection measurements across the Insubric Line from the Sesia Zone into the IVZ near the planned DIVE drilling location at Balmuccia. The primary objective of this seismic survey is to clarify whether the sub-vertical structures prevailing at the surface can be expected to be continuous throughout the uppermost crust, that is, to a depth of ~1 km, or whether there is evidence of major structural flattening or fault-related offsets.
The sub-vertical surficial structures and the high seismic velocities associated with the crystalline rocks present in DIVE's target zone differ strikingly from the layered sedimentary environments, for which the seismic reflection method was originally designed. While seismic reflection profiling has demonstrated significant potential in complex crystalline environments, the acquisition and processing need to be tuned to the corresponding targets, which tend to be non-specular, small-scale, and complex (e.g., Eaton et al., 2003;Malehmir et al., 2012). Correspondingly, the interpretation of the resulting unconventional seismic images is notoriously challenging and tends to lend itself to modelling-based approaches (e.g., Bongajum et al., 2012). In view of this, we have made efforts to acquire seismic reflection data with high spatial and temporal resolution as well as high redundancy, and we geared the processing flow towards enhancing weak seismic events backscattered from complex, non-specular heterogeneities associated with the prevailing sub-vertical structural grain. For the interpretation, we complement the conventional qualitative visual assessment of crustal seismic reflection images with numerical simulations as well as with attribute and geostatistical analyses.

Database
Steeply dipping structures are inherently difficult to image with surface-based seismic reflection surveys, as specular interfaces with dips in excess of ~70 degrees generally do not return any seismic energy back to the surface (e.g., Juhlin et al., 2000). A possibility to address this problem is to focus on the energy backscattered from small-scale, nonspecular heterogeneities associated with larger-scale structures (e.g., Schmelzbach et al., 2008;Tertyshnikov et al., 2015;Khoshnavaz et al., 2016;Schwarz and Krawczyk, 2020). Correspondingly, our acquisition and processing strategies are geared towards collecting high-fold seismic Fig. 1. Geophysical and geological setting of the study area. a) Bouguer anomaly map of the IVZ highlighting the large positive gravity anomaly due to dense rocks at and near the surface. Adapted from Scarponi et al. (2020). b) Lithotectonic map of the central IVZ, framed in white in a). Adapted from Petri et al. (2019). The drilling rig denotes the location of a planned borehole. c) Schematic geological cross-section along the yellow dashed line in b) based on maps and profiles of Hunziker and Zingg (1980), Quick et al. (2003), Berger et al. (2012), and Petri et al. (2019). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) profiles with high spatial and temporal resolution, such that weak backscattered seismic energy can be effectively recorded and enhanced through targeted signal processing and wavefield separation techniques (e.g., Salisbury et al., 2000;Eaton et al., 2003;Schijns et al., 2009;Malehmir and Juhlin, 2010;Juhlin et al., 2010).
The three seismic lines presented in this study, denoted as L1, L2, and L3, were collected along the Boccioleto-Balmuccia and Balmuccia-Isola roads (Fig. 2a). The acquisition and processing methodologies for all three lines are largely identical. In the following, we focus on the data acquired along L1, which crosses from the Sesia Zone in the west into the Insubric Line in the east, because it exemplarily illustrates our methodological approach, quasi-perpendicularly crosses some of the most pertinent structures, and has a significantly higher signal-to-noise ratio (S/N) than L2 and L3.

Acquisition
The data were acquired using a 400 kg weight-drop source (Fig. 2b) in conjunction with 15 Hz vertical-component geophones and a GEODE distributed seismograph system. The acquisition parameters are listed in Fig. 2c. The receiver geometry consisted of a geophone spread with 144 live channels with 5 m station spacing. The receiver spread was rolled along every 24 stations, except for L1, which consisted of a single 48-station roll-on that resulted in a two-spread geometry as shown in Fig. 3a. The source stations were collocated with receiver stations, that is, every 5 m for L1 and L2, which resulted in a dense array of common mid-points or CMPs (Fig. 3b) with a nominal spacing of 2.5 m and a fold of 77 for the 144 live receivers. This fold was increased in the centre of L1 by having overlapping receiver spreads (Fig. 3c). For L3, the source increment had to be increased to 10 m due to time limitations, which corresponds to a nominal fold of 36. The Boccioleto-Balmuccia road section of L1 is only mildly crooked, such that all CMP locations fall within a 30-m-wide swath of the central trend line. Geophones were planted alongside the road and over embankments. Where outcropping solid rock prevented this, 8-mm-diameter holes were drilled into the concrete verge of the road to plant the geophones. Consequently, some receiver station elevations were not vertically collocated with the sources (Fig. 3d), which resulted in erratic first-arrival times (Fig. 4a).
The 400 kg weight-drop source provided signals with a dominant frequency of ~80 Hz (Fig. 4b) and had sufficient energy to reach faroffset receivers, except for the westernmost (~1-30) source positions. In these field records, the first-arriving seismic energy was weak for the easternmost receivers. Lateral variations in signal quality and strength can be seen from west to east in the example field records shown in Fig. 4a, for example, between the dashed vertical white lines. These variations are likely caused by surficial geological variations across the profile and are assumed to be the primary reason for the source energy barely reaching the most distant receivers. This resulted in data with relatively low S/N recorded at these distant channels until the source progressed passed position 30. Additionally, the operation of the weightdrop source was associated with an uncontrollable "bounce", which resulted in ghost shots with onsets at ~300-400 ms. However, due to the high seismic velocities of the study area, these ghost shots, for the most part, do not affect the travel time interval of interest ranging from 0 to ~300 ms. A prominent example of such a ghost shot is shown in Fig. 4a.

Processing
Topographic elevation variations and low-velocity alluvial and/or weathering layers close to the surface are known to significantly distort reflections in hardrock seismic surveys (e.g., Palmer and Jones, 2005;Urosevic and Juhlin, 2007). To alleviate this problem, static corrections were performed prior to pre-stack processing. Firstly, receiver elevation statics were applied to reduce the receivers to a common source-receiver datum corresponding to the Boccioleto-Balmuccia road. Following this, subsurface velocity variations were evaluated using refracted ray traveltime tomography (Fig. 5). The tomographic approach was preferred over standard refraction seismic techniques, such as the plus-minus method (Hagedoorn, 1959) or the generalized reciprocal method (Palmer and Jones, 2005), due to the pronounced lateral velocity variations along the profile, which are clearly observable in the first-arrival travel-time curves (Fig. 5a). Static correction times at each receiver  Horstmann (1981) and Quick et al. (2003). The drilling rig denotes the location of a planned borehole. b) 400 kg weight-drop source used to acquire the seismic data. c) List of key acquisition parameters. location were computed using the inferred surficial velocity structure ( Fig. 5b) down to a depth of 20 m ( Fig. 5c and d), and then applied to the data ( Fig. 5e and f).
Following these static corrections, a targeted seismic reflection processing flow was applied, which is summarized in Fig. 6 and the most essential parts of which are outlined in the following. Conditioning of the pre-stack data involved band-pass filtering (10-20-85-100 Hz) to subdue the airwave, which has a dominant frequency of ~100 Hz, and the subsequent application of spiking deconvolution to balance and enhance the spectral characteristics of data ( Fig. 6a). Spiking deconvolution was successful in boosting the central frequency of the data to ~80 Hz, which corresponds to the upper limit of the previously applied band-pass filter. A non-standard technique was used to remove the direct arrivals and their associated reverberations (Fig. 6b). This was achieved using a 2D alpha-trimmed median filter, which required the flattening of the first-arrivals and the reversal of this procedure after filtering (e.g., Greenwood et al., 2019;Caspari et al., 2020). Please note that this filter also naturally removes the P-wave component of the ghost shots. S-waves, including those associated with the ghost shots, and groundroll were removed using a polygon-type f-k filter. Examples of the pre-stack processing of shot records are shown in Fig. 6c. Here, we see little to no evidence of continuous hyperbolic reflection events that might be associated with shallow gently to moderately dipping structures. Conversely, we see many features that are cross-cutting each other in a chevron-type manner with similar move-outs as the first-arrivals. These opposing dips are analogous to those observed for reflections in VSP data with a vertical receiver line crossing horizontal structures, which represents the transpose to our survey geometry and the prevailing geological structures ( Fig. 6a through 6c). After static corrections, removal of unwanted wavefield components, and CMP-sorting, a corresponding stacked section can be produced. This requires a velocity model to correct for the normal move-out (NMO). NMO velocity analysis is typically based on the assumption of quasihyperbolic move-out associated with a layered subsurface structure. However, in complex, steeply dipping environments, this kind of analysis is not meaningful. Therefore, we visually analysed constant-velocity CMP stacks to determine the stacking velocity, which optimally enhanced the overall strength and coherence of the backscattered seismic energy. While this process is inherently subjective, our key observation was that, while the apexes of the backscattered events changed in time as a function of the stacking velocity, their strength and coherency remained largely constant. Ultimately, we decided on a single, constant stacking velocity of 7000 m/s, which is broadly consistent with evidence from laboratory measurements (e.g., Fountain, 1976). In this context, it is important to note that the use of constant or quasiconstant stacking velocities is quite common in hardrock seismic exploration in response to the small move-outs prevailing in highvelocity crystalline environments (e.g., Eaton et al., 2003).
The resulting stacked section of L1 unsurprisingly does not exhibit any prominent events (Fig. 6d) and, hence, there are also no laterally continuous reflections that could be associated with gently to moderately dipping structures. Conversely, there are many cross-cutting, moderately dipping events, which show coherence over a few traces. These short and discontinuous seismic events, which are quite typical of hardrock environments, are likely to represent seismic energy backscattered from small-scale heterogeneities (e.g., Holliger, 1997;Eaton et al., 2003;Bongajum et al., 2012). There are also weak side reflections coming from ~100 m east of the profile that are back-reflected within the profile at ~750 m (Fig. 6d). These events are interpreted as VSP-type reflections associated with direct waves being reflected from subvertical structures. Finally, imaging and subsequent time-to-depth conversion of the data was performed using Kirchhoff post-stack time migration with a constant velocity of 7000 m/s. In accordance with the presumed backscattered energy present in the stacked section, a small migration aperture of 100 m was used. While the use of a small migration aperture limits the proper imaging of steeply dipping specular reflections, this approach was deemed appropriate due to the absence of such signal characteristics in the stacked data. In turn, this allows to effectively focus the backscattered seismic events and, at the same time, to reduce artefacts, that is, "smiling" effects, in the imaged data. The VSP-type reflections observed in the CMP stack are also visible in the migrated section (Fig. 6e).
Prior to migration there are pervasive, steeply dipping coherent noise events observable throughout the stacked data. The apparent velocity of this noise is ~2000 m/s, thus eliminating the possibility of it being P-waves reflected from vertical boundaries, which would manifest themselves with an apparent velocity of ~3500 m/s. Conversely, the linear move-out of this coherent noise is consistent with an S-wave velocity of ~4000 m/s, which is realistic in the given context. Interpreting this steeply dipping coherent noise as backscattered S-waves and/or groundroll is also consistent with the work of Chopra and Marfurt (2014). Due to their steep dips, these events remain unaffected by migration, which was confirmed by testing migration velocities between 2800 and 7000 m/s specifically for this purpose. To obtain a cleaner final image, we therefore applied a targeted f-k fan filter to the stacked section prior to migration. Extensive tests with varying filter designs demonstrated that this process is robust and largely devoid of artefacts.
Weak, yet pervasive, VSP-type side reflections are still present after migration and cut through the section as shown by the red and yellow lines in Fig. 6e.
The processing flow of L1 outlined above, was also used for seismic lines L2 and L3, and the corresponding stacked and imaged sections are shown in Fig. 7a and b, and Fig. 7c and d, respectively. L2 runs oblique to the geological strike ( Fig. 2) and had an active river gravel quarry operating at its eastern bound. Consequently, the data contained strong and abundant cultural noise, which required very rigorous trace editing. For this reason, the effective average fold of L2 was reduced to ~50 with regard to its original nominal value of 77. The course of L3 is strongly crooked, running sub-parallel to strike (Fig. 2) and on top of Quaternary sediments for the first ~200 m (Fig. 7b). As mentioned above, the shot spacing for L3 had to be increased to 10 m, which reduced its nominal fold to 36. In spite of their lower fold and lower S/N, the overall character of the stacked ( Fig. 7a and b) and migrated ( Fig. 7c and d) sections of L2 and L3 rather closely resembles that of L1 ( Fig. 6d and e).

Interpretation
The resulting stacks and migrated images of all three seismic lines are devoid of laterally coherent reflectors (Figs. 6 and 7), which is consistent with little to no energy being reflected to the surface from larger-scale, specular-type structures. This, in turn, indicates that the Insubric Line, as well as its associated secondary structures, are likely to remain subvertical within the subsurface region imaged by our seismic survey. This first-order assessment is consistent with the fact that the migrated and depth-converted images of L1, L2 and L3, are characterized by local bright lineaments, which subtly change in dip and seem to align sub- vertically. As such, these events could represent the apexes of seismic energy backscattered from local heterogeneities associated with the prevailing larger-scale sub-vertical structural grain (Horstmann, 1981).
In the following, we assess the above hypothesis through the application of attribute and geostatistical analyses of the seismic images. The validity of these approaches, which, as of yet, have not been used in a related context, is tested on synthetic seismic reflection surveys for canonical models of crustal heterogeneity with a sub-vertical structural Fig. 7. CMP stacks of a) L2 and b) L3. Kirchhoff migration images after time-to-depth conversion of c) L2 and d) L3. Please note that the course of L3 is strongly crooked with the first ~200 m running sub-parallel to strike. Fig. 8. a) Gaussian-and b) binary-distributed density models of crustal heterogeneity (Holliger et al., 1993;Holliger, 1996Holliger, , 1997. The density fluctuations of these models are characterized by a von Karman autocovariance function with vertical and horizontal correlation lengths, a z and a x , of 800 m and 200 m, respectively, and a ν-value of 0.3. The P-and S-wave velocities are constant at 7000 m/s and 4200 m/s, respectively. grain.

Synthetic seismic data
To explore the seismic expressions of a sub-vertical structural grain associated with local heterogeneities, we evaluate and analyse synthetic seismic data for canonical crustal models based on Holliger et al.'s (1993) geostatistical conceptualization of Ivrea-type lower crust (Fig. 8). We consider both a standard Gaussian-distributed stochastic model as well as a binarized version thereof, which allows us to explore the endmember-type seismic responses of smoothly varying and abruptly changing material properties along sub-vertical structures (Fig. 8). The density fluctuations in both models are characterized by a so-called von Karman autocovariance function (e.g., Tronicke and Holliger, 2005) with vertical and horizontal correlation lengths a z and a x of 800 m and 200 m, respectively, a Hurst number ν of 0.3, and a constant mean density of 2800 kg/m 3 . Both the P-and S-wave velocities are kept constant at 7000 m/s and 4200 m/s, respectively, to facilitate the subsequent imaging and to reduce associated biases. The Gaussiandistributed model is characterized by standard deviation of the density fluctuations of ~100 kg/m 3 , while the density values in the binarydistributed model are either 2600 kg/m 3 or 2800 kg/m 3 . As such, the average impedance contrasts in the binary-distributed model are significantly higher than those of its Gaussian-distributed counterpart. The considered autocovariance model as well as its correlation lengths and ν-value are consistent with the results of Holliger et al. (1993) inferred from the combined stochastic analysis of geological maps and rock physical properties for some key locations in the IVZ. The binarydistributed model can be regarded as a first-order approximation of a complex interlayering of mafic and felsic rocks locally prevailing in the IVZ (Holliger et al., 1993), while the Gaussian-distributed model emulates compositionally less heterogeneous upper crustal granitic and gneissic rocks (e.g., Holliger, 1996Holliger, , 1997 and, hence, can be expected to rather representative of the gneisses and mylonites prevailing in the Sesia Zone and the Insubric Line, respectively.
For these models, we then generated synthetic seismic reflection surveys using a staggered-grid finite-difference solution of the elastic wave equation, which is fourth-order accurate in space and secondorder accurate in time (e.g., Levander, 1988). All key acquisition and processing parameters were chosen to emulate those of the field data. The source is a Ricker wavelet with a centre frequency of 80 Hz. Figs. 9 and 10 show the resulting stacked and migrated synthetic seismic sections. The data for the binary-distributed model exemplarily illustrate the potential of the chosen acquisition and processing strategy for detecting and imaging seismic energy backscattered from local heterogeneities associated with larger-scale sub-vertical structures characterized by strong impedance contrasts. This backscattered energy is also present, albeit in a much less prominent manner, in the synthetic data for the Gaussian-distributed stochastic model. In this context, it is important to note that, in terms of the overall characteristics, both the stacked and migrated field data of L1 bear a conceptual similarity with their synthetic counterparts for the Gaussian-distributed model. This observation is consistent with the fact that L1 is primarily located in gneissic and mylonitic rocks, for which the fluctuations in the material properties are indeed expected to be continuous and quasi-Gaussian (e. g., Holliger, 1996Holliger, , 1997.

Attribute analysis and geostatistical inversion
To enhance the seismic energy backscattered from heterogeneities associated with sub-vertical structures, we subjected the migrated versions of both synthetic and field data to energy-based attributed analysis (Fig. 11). This attribute is commonly used in sedimentary environments to assist the detection and characterization of lateral structural discontinuities associated with, for example, faults and chimneys (e.g., Avseth et al., 2010). Here, we explored the potential utility of this attribute in complex hardrock environments for revealing the coherence between individual "bright lineaments" associated with seismic energy backscattered seismic from local heterogeneities associated with the overall structure grain. The backscattered energy was calculated using a moving window with a length of 60 ms along each trace and defined as the sum of the squared sample values within the window normalized by the number of samples. Fig. 11 shows a comparison of the energy images of synthetic seismic data for Gaussian-and binary-distributed models with that for the field data recorded along L1, L2, and L3. This illustrates the consistency amongst the observed data as well as their similarity in character with the synthetic data in general and those for the Gaussiandistributed model in particular, which, in turn, qualitatively supports the hypothesis outlined above and, thus, the prevalence of the subvertical structures throughout the imaged part of the upper crust.
To further corroborate and quantify this interpretation, we proceeded to analyse the migrated synthetic and field data using the method of Irving et al. (2009), which relates the geostatistical properties of images of the backscattered wavefields to those of the scattering media. Specifically, this method allows to estimate the structural aspect ratio of the underlying heterogeneity and, to a lesser extent, also its complexity, which for the von Karman autocovariance model is quantified by the Hurst number ν. While this technique has been applied successfully to heterogeneous sub-horizontally structured environments (e.g., Irving et al., 2009;Irving and Holliger, 2010;Scholer et al., 2010;Xu et al., 2020), its applicability in the presence of sub-vertical structures was as of yet unproven. In the following, we therefore first assess Fig. 9. CMP-stacked synthetic seismic reflection data obtained for the a) Gaussian-and b) binary-distributed heterogeneous crustal models shown in Fig. 8. Y. Liu et al. the viability of this approach for our synthetic data and subsequently apply it to the observed data. Fig. 12 shows the corresponding validation for the migrated synthetic seismic images for the Gaussian-and binary-distributed models as well as the results obtained for the field data recorded along L1, L2, and L3 based on a Monte-Carlo-type inversion approach (e.g., Xu et al., 2020). Please note that the first ~200 m of L3, which run sub-parallel to strike, have been excluded from this analysis. We see that the resulting statistics for the synthetic data ( Fig. 12a and b) are largely consistent with the ratio of the vertical to horizontal correlation lengths a z /a x of 4 for the underlying models. The results for the Hurst number ν are, as expected from previous studies (e.g., Xu et al., 2020), less well constrained, but nevertheless in the correct range. These results thus seem to confirm the validity of the method of Irving et al. (2009) even in the presence sub-vertical structures. The corresponding analyses of the field data also point to a structural aspect ratio a z /a x that is significantly larger than unity as well as to a rather small ν-value. While the former agrees with the presumed dominance of sub-vertically aligned heterogeneities, the latter is consistent with the observation that the stochastic distribution of elastic properties in crystalline rocks seems to be universally characterized by low ν-values (e.g., Holliger, 1996Holliger, , 1997). An interesting observation is that the character of the histogram of the a z / a x -values for L2 differs from that for the synthetic data as well as that for L1 and L3, which may be due to the fact the L2 runs oblique to strike and has a lower S/N than the two other lines.

Dip angle analysis
The above qualitative and quantitative analyses indicate that the outcropping sub-vertical structures are likely to persist throughout the upper crustal section imaged by our seismic data. The question that still needs to be addressed concerns the average dip of the structural grain, which, provides information as to whether the dominant dips observed at the surface remain more or less constant or undergo significant Fig. 10. Migrated and depth-converted images of the stacked synthetic seismic reflection data obtained from the a) Gaussian-and b) binary-distributed heterogeneous crustal models shown in Fig. 8.   Fig. 11. Energy-based attribute analysis of the migrated depth images of the synthetic data for a) Gaussian-and b) binary-distributed heterogeneous crustal models as well as for the field data recorded along c) L1, d) L2, and e) L3. Fig. 12. Results of the Monte-Carlo-type geostatistical inversion based on the method of Irving et al. (2009) of the migrated depth images of for vertical-to-horizontal aspect ratio (left column) and the ν-value (right column) for the synthetic seismic data based on a) Gaussian-and b) binarydistributed heterogeneous crustal models and the field data recorded along c) L1, d) L2, and e) L3. The mean values for a z /a x in a), b), c), d, and e) are 4.5, 4.9, and 3.9, 4.7, and 3.1 respectively. The mean values for ν in a), b), c), d), and e) are 0. 25, 0.25, 0.26, 0.32, and 0.32, respectively. changes at depth. We found that this problem can be approached by analysing the secondary lobes of the autocorrelations of the energy images (Fig. 11), which reveal the larger-scale alignment of the backscattered energy across the entire section. To assess and validate this novel approach, we again generated synthetic seismic reflection data for our stochastic models with the same structures as before (Fig. 8), but inclined by 10, 20, and 30 degrees with regard to the vertical. Figs. 13 and 14 show the corresponding Gaussian-and binary-distributed seismic models together with the 2D autocorrelations of the corresponding energy images of the synthetic seismic data. While there is some interpretational leeway, the trend of secondary lobes of the autocorrelations is clearly indicative of the overall structural dip. Fig. 15 shows the autocorrelation of the energy image of field data collected along L1, L2, and L3. Note that the first ~200 m of L3, which run sub-parallel to strike (Fig. 3), have again not been considered in this analysis. The autocorrelations of the energy images of L1 and L3 exhibit prominent side lobes whose trends point to steep dips of the order of 10 degrees in the eastwest direction. This is consistent with the predominant dips prevailing at the surface (Horstmann, 1981) and seems to confirm that the outcropping structural grain seems to persist throughout the seismically imaged part of the upper crust. Conversely, the result for L2, which runs oblique to strike over its entire length and has a lower S/N than L1 and L3, is more ambiguous due to the absence of well-defined side lobes.
Finally, Fig. 16 provides a synoptic consolidation of the results of this study in form of an overlay-type representation of the seismic depth images and the associated energy-based attributed of L1, L2, and L3, juxtaposed with the local surface geology mapped by Horstmann (1981) and Quick et al. (2003). Also shown are the dips of the predominant structures, namely the mylonitic border zone of Insubric Line crossed by L1 and the gabbro/peridotite contact near Balmuccia crossed by L3. For both of these sub-vertical structures, there is a clear correlation with changes in seismic character whose vertically continuous and laterally discontinuous nature is highlighted by the energy-based attributes. This representation of the seismic data is consistent with our previous qualitative and quantitative analyses, which indicate that the sub-vertical structural grain prevailing at the surface is likely to continue at least throughout the first one to two kilometres of subsurface imaged in our survey.

Discussion and conclusions
We have presented the acquisition, processing, and interpretation of a high-resolution seismic reflection survey crossing the Insubric Line from the Sesia Zone into the IVZ. The objective of this survey, which consists of three ~1-km-long lines of variable S/N, was to characterize the first ~1 km of the subsurface in preparation for an ICDP drilling project aiming to explore lower continental crust exposed in the IVZ. Given that the structures exposed at the surface are sub-vertical, the key challenge for this project was to reliably assess whether the outcropping structural grain is continuous at depth or whether there are any indications for major flattening or shallow-angle cut-offs. Given this scenario, specular-type reflections were unlikely and the backscattered energy, which seemingly manifests itself as local bright lineaments in the seismic images, was tentatively interpreted as originating from small-scale heterogeneities associated with the sub-vertical large-scale structural grain.
To test this hypothesis and to corroborate the associated qualitative interpretation of the seismic images, we have generated synthetic seismic reflection surveys for sub-vertically structured canonical models of crustal heterogeneity. The acquisition and processing parameters of these synthetic seismic surveys emulate those of the field data. We consider both Gaussian-and binary-distributed heterogeneities for our models. While the stacked and imaged synthetic data for the Gaussiandistributed models bear a clear resemblance with the field data, the synthetic data for binary-distributed models exemplarily illustrate that the recorded backscattered energy originates from local heterogeneities Fig. 13. a) Gaussian-distributed crustal models with the same structures as in Fig. 7 but dipping with 10, 20, and 30 degrees with regard to the vertical. b) 2D autocorrelations of the energy images of corresponding synthetic seismic data. Y. Liu et al. associated with the overall sub-vertical large-scale structural grain. This interpretation was corroborated by performing energy-based attribute analysis as well as by a geostatistical inversion of the depth images of the observed and synthetic seismic data with regard to the structural aspect ratios of the backscattering structures. For the latter, we first needed to demonstrate the validity of the underlying method, which was originally conceived for sub-horizontal structured environments, in the presence of sub-vertical structures. Although, this geostatistical analysis is less sensitive to the complexity of the underlying medium than to its structural aspect ratio, the corresponding results are consistent with the ubiquitous and seemingly universal characteristics of crystalline rocks. Finally, we illustrate that the alignment of the side lobes of the 2D autocorrelation of the seismic energy images allows to assess the overall dip of the prevailing larger-scale structural grain, which for the observed data is consistent with the sub-vertical east-west dip of the predominant structures mapped at the surface.
The qualitative and quantitative analyses of our seismic reflection survey thus indicate that the sub-vertical structural grain exposed at the surface is likely to prevail throughout the imaged part of the upper crust. Beyond addressing this specific objective, the multi-faceted and unconventional approach pursued in this study provides interesting insights and potentially new perspectives with regard to the acquisition, processing, and interpretation of high-resolution seismic reflection data in crystalline terranes and their capacity for providing information on complex structures whose dips are too steep for deterministic imaging. Conversely, it is important to bear in mind that this seismic reflection survey violated the method's founding assumptions and, while the results of our unconventional interpretation approaches seem interesting Fig. 14. a) Binary-distributed crustal models with the same structures as in Fig. 7 but dipping with 10, 20, and 30 degrees with regard to the vertical. b) 2D autocorrelations of the energy images of corresponding synthetic seismic data. Fig. 15. 2D autocorrelation of the energy images of the field data recorded along a) L1, b) L2, and c) L3. Please note that the westernmost ~200 m of L3, which run sub-parallel to strike, have not been included in the analysis. and promising, we are unable to authoritatively validate and generalize them at this point. A critical question in this context, and indeed for most, if not all, 2D seismic reflection profiles, concerns the importance of 3D effects. This is particularly pertinent in our case, as, although the strike and dip of the larger-scale geological structures are quite stable along the survey, this quasi-2D assumption cannot be extended to the smaller-scale heterogeneities, which are superimposed on these larger structures, and which represent the very target of our interpretation efforts. Rigorously addressing the associated question regarding the importance of 3D scattering effects on our results would require extensive 3D finite-difference-type elastic modelling for a wide range of scenarios. While this is clearly beyond the scope of the current work, it represents an interesting and pertinent research problem in its own right, which ultimately will need to be addressed in order to assess the reliability and robustness of the proposed interpretation approaches. At present, our admittedly quite speculative hypothesis in this regard is that, due to their inherent 3D nature, out-of-plane scattering effects would have a tendency of reducing, rather than enhancing, the subvertical alignment of backscattered seismic events with regard to the considered 2D reference case. This, in turn, would imply that we tend to underestimate, rather than overestimate, the ratio of the vertical to horizontal correlation lengths a z /a x.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 16. Presentation of the seismic data of L1, L2, and L3 relative to the geologic mapping of Horstmann (1981) and Quick et al. (2003). The seismic data are displayed in greyscale overlain with the colorscale energy-based attributes from Fig. 11. Please note that the first ~200 m of L3, which run subparallel to strike are not shown. The major structural dips mapped at the surface (grey circles) are denoted within the sections by solid red lines, while their presumed continuation at depth is denoted by dashed red lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)