Coherent diffraction imaging for enhanced fault and fracture network characterization

Faults and fractures represent unique features of the solid Earth and are especially pervasive in the shallow crust. Aside from directly relating to crustal dynamics and the systematic assessment of associated risk, fault and fracture networks enable the efficient migration of fluids and, therefore, have a direct impact on concrete topics relevant to society, including climate-change mitigating measures like CO2 sequestration or geothermal exploration and production. Due to their small-scale complexity, fault zones and fracture networks are typically poorly resolved and their presence can often only be inferred indi5 rectly in seismic and ground-penetrating radar (GPR) subsurface reconstructions. We suggest a largely data-driven framework for the direct imaging of these features by making use of the faint and still often under-explored diffracted portion of the wavefield. Finding inspiration in the fields of optics and visual perception, we introduce two different conceptual pathways for coherent diffraction imaging and discuss respective advantages and disadvantages in different contexts of application. At the heart of both of these strategies lies the assessment of data coherence, for which a range of quantitative measures is introduced. 10 To illustrate the approaches versatility and effectiveness for high-resolution geophysical imaging, several seismic and GPR field data examples are presented, in which the diffracted wavefield sheds new light on crustal features like fluvial channels, erosional surfaces, and intricate fault and fracture networks on land and in the marine environment.

With only few exceptions (e.g. Landa et al., 1987;Heincke et al., 2006;Dell et al., 2019), the potential of quantitative coherence measures for directly forming noise-robust, contrast-rich images remains largely unexplored. Building on recent 60 advances in adaptive processing and weak-wavefield enhancement, we present a strategy for seismic and GPR diffraction imaging that makes direct use of wavefield coherence for scatterer detection. After a brief elaboration on typical characteristics and unique properties of diffraction phenomena, we introduce two different means of reconstructing a scatterer with the help of coherence measurements. Underpinning both these pathways we introduce generalized coherence measures and systematically investigate their tolerance with respect to imperfect, i.e. noisy, sparse, or incomplete data and make suggestions with respect to 65 their applicability. Concluding community-spanning seismic and electromagnetic examples suggest that coherent diffraction imaging not only leads to overall highly-resolved subsurface reconstructions, but also directly and reliably highlights smallscale erosional features, faults and fractures.

Wave diffraction
Diffraction can loosely be defined as a waves ability to enter shadow zones, which are forbidden regions in geometrical optics. 70 More precisely, diffraction occurs when a wavefield encounters a relevant property change that has a local curvature of or below the wave length (Born and Wolf, 2013). Thus, diffraction is a scale-spanning phenomenon that is only predicted and fully captured in a wave theoretical framework. To provide some intuition, Figure 1 illustrates some of the key properties of diffractions and how they can be of use for geophysical subsurface imaging (for more details, see Schwarz, 2019a). As arguably the first rigorous experimental evidence, Young's slit experiments concluded that when light hits a small enough opening in a 75 screen, an intricate interference pattern appears on a second screen. The geophysical analogue of such an experiment is shown in Figure 1(a), where a small-scale heterogeneous object or a gap at an interface acts as a secondary source. The acquisition surface at zero depth can be viewed as a screen, where the wavefields are captured by seismometers or electromagnetic antennas. The mere occurrence of such a secondary wavefield (for clarity the primary field is not displayed) is already indicative of the presence of a small-scale structural change underneath our acquisition surface, which is why it is frequently suggested, that 80 the detection and localization of such a structure potentially implies super-resolution imaging, i.e. the inference of structural features of spatial extent beyond the Rayleigh limit (Khaidukov et al., 2004).
Because their wavefronts are principally indistinguishable from an actual source located at the scatterer location, seismic diffractions share striking similarities with passive sources (Li et al., 2020). In Figure 1(b) the surface projections of the diffracted wavefield shown in (a) are displayed as a function of time. Generally, although this is not precisely true in realistic 85 media, diffractions have approximately hyperbolic shape. In contrast to reflections, diffractions are not constrained by Snell's law and, thus, radiate uniformly in all directions. As a result, diffracted wavefields provide improved illumination and encode highly resolved structural information, but also rapidly decay with increasing distance from the scatterer. Closely related to diffraction is the concept of interference, which likewise is a pure wave phenomenon. Interference is mentioned here for two reasons. First, it explains the transitional regime and provides a notion of resolvability (Figure 1(c)), and second, it helps 90 to appreciate the need and possibility for wavefield separation, in particular, when highly reflective media like sedimentary basins are considered (Figure 1(d)). When a sufficiently large amount of scatterers is present, individual diffractions become hardly distinguishable. Essentially all diffraction separation strategies rely on dense spatial sampling at the surface. The highresolution and high-illumination component of diffractions, which bear unique imaging potential, can only be unlocked, if spatial aliasing can be prevented (Schwarz, 2019a).

95
Aside from illustrating lateral resolvability, Figure 1(c) also visualizes the underlying principle of Kirchhoff migration (Schneider, 1978). As will be more thoroughly explained in the following section, Kirchhoff migration is a manifestation of the Huygens-Fresnel principle, which states that any arbitrary wavefield can be thought of as being composed of infinitely many elementary wavefields. The envelope of these elementary, locally excited waves then forms the transmitted or reflected arrival.
Diffractions can be viewed as physically resolved manifestations that are picked out by small-scale disturbances, as e.g. caused 100 by faults. Figure 1(d) shows a small excerpt of a subsurface model that mimics the Sigsbee escarpment in the Gulf of Mexico.
Aside from structural features related to the top of salt or a fault cutting through the sedimentary strata, wave diffraction is also caused by stair stepping in the discretised model used for the finite-difference simulations. Despite their unintentional introduction, the resulting pervasive diffracted wavefields nicely illustrate the transition from diffraction to reflection and why a successful separation of these contributions remains a challenge to confront. High-resolution imaging aims at back-tracing 105 these weak diffracted contributions to their origin.

Coherent diffraction imaging
In conventional diffraction imaging, wavefield separation is either performed before or during migration (see, e.g., Fomel et al., 2007;Moser and Howard, 2008). Independent of what type of migration algorithm is used, the result commonly comprises a wavefield image that contains amplitude and phase information. While the preservation of phase information in the reconstruc-110 tions is principally desirable, there exist several shortcomings of conventional migration schemes, in particular when imperfect data and weak signals such as diffractions are concerned. In coherent diffraction imaging we seek to directly incorporate wavefield coherence in the imaging workflow to help overcome these limitations. It is generally interesting to note that when optical images are concerned, we only perceive intensities, and wavelength information is encoded in the colours of the visible light, to which the eye is sensitive. Following this intuition, we argue that coherence measures, to some degree, mimic intensities and, 115 therefore, seem principally suited for the construction of structural images. As was briefly illustrated in the previous section, diffractions are coherent and, just like reflections, can benefit from coherence arguments.
The imaging problem can generally be divided into two domains -the data and the image. Migration-type imaging, just like an optical lens, seeks to directly utilize the former to arrive at the latter. In both, data and image space, coherence arguments can be invoked ( Figure 2). While the systematic assessment of coherence in data space has a long and successful history, using 120 coherence measures, with some few exceptions, have not been utilized to their full potential, when the image space is concerned.
With intuition from the field of optics, in order to properly differentiate between these two philosophies, we refer to starting in data space, i.e. with the observations, as projection, whereas the image-centric approach will be denoted by focusing. Both mindsets have in common that we use the data to construct an image and both are amenable to improvements when some form of coherence measure is introduced. Consequently, coherent focusing evaluates wavefield coherence during the gathering stage 125 in image space, whereby coherent projection first evaluates data coherence and then back-projects with help of the extracted information. Wave equation migration (focusing) and time-reversal imaging (projection) can be viewed as the most capable end members of these two branches which find themselves powerfully combined in reverse-time migration for reflection imaging (Baysal et al., 1983). While fully honoured wave propagation physics becomes important in sufficiently complex scenarios, the unique flexibility of Kirchhoff migration and its intimate relation to wave diffraction provides unique opportunities for the 130 imaging of scatterers (Moser and Howard, 2008). Although wave propagation is abstracted by high-frequency approximations of limited validity, the use of only kinematic information lets the developed framework be readily applicable for both, seismic and GPR measurement campaigns. tion. While projection-type imaging schemes (Section 3.3) start directly in data space, focusing techniques (described in Section 3.2) typically are image-centred. Although diffraction separation prior to imaging forms a central ingredient, emphasis is given to the reconstruction steps.

Measuring coherence
As illustrated in Figure 1, coherence is an observable contained in densely acquired waveform data. It can be assessed in depth 135 (Figure 1(a)) and time (Figure 1(b)), i.e. in image and in data space. Coherence can be viewed as a set of correlations that are connected by the temporal or spatial delays arising from the shape of the propagating wavefront. If data are not acquired densely at the surface, these group correlations cannot be reliably tracked and connected any more. A well-appreciated way of assessing wavefield coherence is to perform directional data summations within a predefined time window. If we denote the spatio-temporal waveform data recorded at the surface (at position x and time t) by D(x, t), we can write the summed amplitude at point ξ, estimated within a confined aperture spanned by all x, as (1) where t(ξ, x) is the traveltime surface corresponding to the wavefront. If this traveltime surface describes an actual event (compare Figure 1(b)), the summation result B(ξ) shows increased amplitudes, while for uncorrelated noise or a wrong or inaccurate choice of t(ξ, x), the amplitudes are smaller. Equation (1) can be used to systematically suppress uncorrelated noise 145 or undesired interfering coherent energy, but does not lend itself well for an automated analysis of wavefield coherence. In addition, the summed wavefield, like the data itself, encodes phase information resulting in positive and negative values, which complicates interpretation. In analogy to optics, we will refer to expression (1) as the beam amplitude or beam, which follows the physically intuitive convention in earthquake seismology (Rost and Thomas, 2002). To arrive at a more robust quantity that can act as a cost function in an optimization scheme, the beam energy can be approximated as follows where τ is a small time window in which vertical summation is performed. As a rule of thumb, it should have approximately the size of the considered signal's wavelength. The beam energy E(ξ) takes only positive values, but does not precisely correspond to the beam's energy but rather is proportional to it. In earthquake seismology, equation (2) is investigated routinely in slowness and back-azimuth analysis (Rost and Thomas, 2002). If we consider the total energy contained in the investigated portion of 155 the wavefield, a similar proportionality holds Expression (3) lets us arrive at an upper bound, as coherent summation has to honour energy conservation. The wavefields thus is a normalized quantity and was demonstrated to be an ideal candidate for the automation of coherence analysis (Taner and Koehler, 1969;Neidell and Taner, 1971). The quantity N x indicates the total number of traces at all recording locations x that fall within the considered aperture. For perfectly coherent signal S approaches 1, whereas for fully uncorrelated noise it takes values close to 0. If data are very noise contaminated or other contributions interfere strongly with the event investigated, it can be useful to abstract the waveform data before processing (Li et al., 2020). One such abstraction constitutes the polarity-165 honouring n-th root of the signal One main advantage of the n-th root abstraction is that in contrast to other means like kurtosis, the transformed data still retain their polarity, which allows for destructive interference to occur. Insertion of expression (5) for D(x, t) in equations (1)-(4) then leads to n-th root versions of the beam amplitude B n (ξ), the beam energy E n (ξ), and semblance S n (ξ) (compare Schwarz, 2019a). For n = 1 all of these expressions reduce to their conventional analogues, which is why in the following, we will refer to different versions of each quantity by their order n. A systematic investigation of the above coherence measures will be carried out in the following two subsections. Nevertheless, it can already be stated that the beam energy (2) bears a strong resemblance to a wavefields intensity, which likewise is sensitive to the absolute amplitude of a signal. As a consequence, it ascribes higher values to stronger, more energetic contributions (stronger scatterers appear brighter). Conversely, semblance 175 represents an energy ratio and, owing to its normalization, coherence is detected independent of signal strength.

Imaging by focusing
As indicated earlier, imaging by focusing can conveniently be based on Kirchhoff's diffraction integral, which in practice, like equation (1) reduces to a discrete sum. In a more physical sense focusing can be viewed as a special form of constructive interference, in which different measurements of the same coherent wavefield are superposed at the image location ξ. This 180 image point represents either a sample in the spatio-temporal focussed image ξ = (x 0 , t 0 ) corresponding to time migration, or fully spatial reconstruction with a depth axis ξ = (x 0 , z 0 ) (Schneider, 1978;Etgen et al., 2009). As we are primarily interested in robust structural images of small-scale crustal features, here, in line with other authors (Khaidukov et al., 2004;Moser and Howard, 2008), we neglect individual amplitude weighting which is normally applied in true-amplitude migration. As a consequence, the discrete version of the Kirchhoff integral is equivalent to the beam amplitude (1). Although all the presented 185 findings likewise translate into depth, for the sake of simplicity and in order to stay consistent with the field data examples for which detailed velocity information was scarce, in the following, we will only consider spatio-temporal focusing.
In order to arrive at a full equivalence of the un-weighted Kirchhoff integral with the beam amplitude, for every image point ξ = (x 0 , t 0 ), the traveltime response of a diffraction needs to be inserted into equation (1). While for depth migration the accurate traveltimes of a diffracted wavefront are computed by means of ray tracing or eikonal solvers, time imaging typically 190 relies on analytical, closed-form approximations whose validity is restricted to mild levels of lateral heterogeneity. A popular assumption is that traveltime moveout is approximately hyperbolic, which corresponds to circular wavefronts in an effective replacement medium (e.g. Schwarz and Gajewski, 2017b). This effective medium is described by the root-mean-square velocity v rms , which generally is a function of the image point. In compact notation, the multi-channel diffraction traveltime t diff can be written as a sum of a source and a receiver leg connected by the shared image point Figure 3 illustrates with a synthetic example, how the previously discussed coherence measures compare for high-quality, i.e. densely and regularly sampled data with a low noise level, severe noise contamination, data sparsity resulting in spatial aliasing, and incomplete (single-sided) observations. Without loss of generality, alongside the un-weighted Kirchhoff reconstruction, i.e. the beam amplitude (1) with t = t diff , only the first and the 10-th order of the beam energy (2) and 200 the semblance coefficient (4) is displayed. For the high-quality data (top row in Figure 3) all considered coherence measures arrive at an accurate reconstruction of the two scatterers, while in the case of data insufficiencies especially the conventional Kirchhoff-type reconstruction with the beam amplitude suffers from strong imaging artefacts. In Kirchhoff migration, diffractions are naturally favoured in that in contrast to reflections, the summation trajectory is not only tangential, but fully coincides with the event. Because the two diffractors are located in different depths and at different lateral positions, they are not equally 205 well imaged in all scenarios. Especially when only the left-sided incomplete observation is concerned, the right diffractor remains poorly resolved with the beam amplitude, the conventional beam energy, and the conventional semblance norm.
As a typical challenge in diffraction imaging constitutes the presence of residual reflected energy after separation, we included the response of a planar reflector in the deeper part of the model. While this strong but undesired contribution is only mildly suppressed in the conventional Kirchhoff reconstruction, the coherence images are mostly reflection-free with only 210 weak residual energy remaining for the conventional beam energy and the semblance norm. Both, the 10-th root version of the beam energy and the 10-th root semblance lead to diffraction focusing results of consistently high quality in all evaluated scenarios. Following from this systematic analysis and because its normalized character that favours weak and incompletely sampled wavefields, we conclude that the n-th root semblance can be expected to be the most robust candidate for diffraction imaging with imperfect data. The strong suppression of reflected energy furthermore suggests that even without competitive 215 diffraction separation first diffraction-enhanced images might be gained in moderately reflective media. Taking the n-th root of the amplitude as suggested in expression (5) has the effect of making coherent arrivals of different strength more comparable.
While the suggested value of 10 results from experience with a variety of data configurations, it can be shown that this equalization in amplitude typically saturates for a reasonably low n already. In principle, the problem of finding a suitable root order could be phrased as an optimization problem, driven by the amplitude content of the data. However, a fixed value of 10 was 220 shown to be successful in bridging several orders of magnitude and, therefore, is deemed a reasonable choice in most practical scenarios.
Very similar to the anisotropic radiation characteristics of passive seismic sources, for diffraction off edges and structural steps, as they likely also occur at fault zones, the polarity of diffractions can change at the stationary point, which would lead to a bi-modal reconstruction. To account for this radiation pattern, the conventional coherent focusing result can be augmented 225 with its phase-reversed version ( Figure 4). Every data point is once treated as a potential stationary point at which an artificial phase reversal is performed before evaluating the coherence measure by means of equations (2) or (4). Both results, the one gained without reversing the phase, and the one for which the phase is reversed, are compared and the higher value contributes to the augmented image. The augmented counterparts of the beam energy and the semblance norm are insensitive to the occurrence of a phase reversal, resulting in a more stable and only slightly less resolved reconstruction. Additionally, it can act 230 as a data-driven migration weight to suppress artefacts in Kirchhoff-type diffraction focusing (compare the rightmost column in Figure 4).

Imaging by projection
An alternative to diffraction focusing constitutes imaging by projection. The main mindset underlying this second strand of diffraction imaging is to start investigations in data space and to use the extracted information to arrive at an image. In coherent 235 diffraction imaging, coherence analysis is carried out in the data domain to locally enhance and physically characterize waveform similarities that can be exploited for imaging. The previously discussed coherence measures can be readily employed.
However, instead of using the reconstruction-centred image point parametrisation, the emerging diffracted wavefront is locally characterized along the entire event, resulting in a transformation of the data volume into coherence and wavefront attribute maps. A data-centred 2D analogue to the double-square-root equation (6) can be expressed in terms of the local tilt angle α 0 240 and curvature radius R 0 of the wavefront as it is observed at location x 0 on the acquisition surface  Figure 3. The augmentation of a phaseconsistent and a phase-reversed version of each considered coherence measure leads to a reconstruction that is slightly less resolved laterally, but proves to be insensitive to a sign change in amplitude for diffraction off asymmetric objects like edges. In addition, the normalized semblance norm can be used as a fully data-derived amplitude weight to suppress artefacts or noise in migrated images. Höcht et al., 1999;Schwarz and Gajewski, 2017b). It is important to note that in contrast to the process of image formation by focusing, here, the discussed coherence arguments are evaluated within a local data aperture and assigned to the central data point within this aperture. So in contrast to equation (6), (x 0 , t 0 ), just like the 245 summation process itself, now also resides in data space. The actual reconstruction then consists of a subsequent mapping of every locally coherent data point to a point in image space. In analogy to focusing-based time migration, analytical mapping equations can be gained, by evaluating the stationary point (x apex , t apex ) corresponding to the apex position of each individual local diffraction traveltime curve. For expression (7), the projection corresponds to a mapping from each data point (x 0 , t 0 ) to where v 0 denotes the locally constant near-surface velocity. Generally, by convention, we refer to near-surface quantities, i.e. quantities that relate directly to the registration, with the subscript 0. Equations (8) in conjunction with expression (7) and one of the positive-valued coherence measures can be used to set up an optimization problem to arrive at an approximate but 255 fully data-driven reconstruction of the subsurface scatterer distribution (Fomel, 2007;Schwarz et al., 2014;Bonomi et al., 2014). In addition to considering a near-surface projection, wavefront slopes and curvatures can also be estimated using the assumption of an effective replacement medium, which, like in equation (6) is defined by the root-mean-square velocity. A so-called osculating equation connecting the near-surface projections and effective medium properties was first established by de Bazelaire (1988) and generalized by Schwarz and Gajewski (2017b). For the sake of simplicity, only the 2D versions of the 260 Figure 5. Results of data-driven coherence analysis (left) and subsequent projection (right) carried out for the same synthetic test case as in In Figure 5 the maximized coherence for every data point as well as the corresponding projection results for the same synthetic test as in Figure 3 and Figure 4 are presented. While all four considered coherence measures perform equally well when the mere detection of coherent energy is concerned, noticeable differences can be observed with regards to overall 265 strength and the handling of interference. Both versions of the beam energy turn out to be not affected by the presence of the strong conflicting reflected event, whereas semblance, due to its intrinsic normalization, reveals to be suited to estimate the coherence of weak and energetic arrivals equally well. In contrast to the results gained with the conventional definition of the semblance (n = 1), the 10-th root version does not suffer from interference effects and, similar to the preceding focusing analysis, performs consistently well for all data points. Although all four displayed coherence measures can be used as a cost 270 function in data-driven coherence analysis, the use of normalized quantities allows for efficient and intuitive thresholding in subsequent processing steps, which, again, lets the 10-th root semblance appear as the most favourable candidate. Similar to coherent focusing, reflections are naturally suppressed by being projected in a diffuse manner.
As will be commented on in more detail in the discussion section, the data-centric nature of projection favours automation and macro-model-independent imaging. However, a major difference between data and image space constitutes the occurrence 275 of intricate interference phenomena, in which multiple wavefields conflict with each other. This can be viewed as an efficient means of compression, but the decoding (i.e. the separation) of interfering wavefield components remains a challenging and computationally demanding task (compare the degradation of the conventional semblance norm in Figure 5). In contrast to that, the image-centric character of coherent focusing naturally avoids these complications. For this reason, without loss of generality, the challenging data examples considered in the following section will all be based on coherent focusing.

Applications
The occurrence of diffraction phenomena is linked to the pre-dominant wavelength. Consequently, just like faults and fractures themselves, diffractions can be observed on essentially all scales probed in geophysical investigations of the upper crust.
As with conventional Kirchhoff migration, there exist natural limitations of the suggested strategies for coherent diffraction imaging. However, with a range of ambitious field data examples we seek to demonstrate that rich diffracted wavefields exist 285 in essentially all datasets and become assessable with the presented robust coherence arguments. Following the quantitative analysis of the different coherence measures discussed in the previous section, the images presented in the following were without exception generated with the 10-th root semblance norm. It may however be noted that the other coherence measures might have led to reconstructions of comparable quality. With exception of the 3D seismic land data example, all diffraction images were formed by augmenting sub-images with and without radiation pattern correction. Despite the depiction of the gen-290 eralized workflow for coherent diffraction imaging in Figure 2, the presented images did not experience any image processing, but relied on a preceding separation step to suppress reflected energy from the considered pre-and poststack data (Schwarz, 2019b). While ground-truth reconstructions were not available, a systematic comparison with Kirchhoff-migrated images is provided for some of the examples.

295
The first data example consists of a marine multi-channel seismic field dataset acquired by TGS in the context of hydrocarbon exploration in the Levantine Basin, Eastern Mediterranean (e.g. Reiche et al., 2014). In contrast to the other data examples, this dataset includes source-receiver offsets of up to more than 7 km, which makes it well suited for conventional migration and inversion. The captured geology off the coast of Israel is primarily dominated by the presence of salt and connected processes in the upper crust (Krijgsman et al., 1999;Gradmann et al., 2005;Netzeband et al., 2006). Most notably, a large, laterally elongated

Single-channel seismic imaging in the Aegean Sea
As a second example we consider a single-channel marine seismic dataset that was acquired near the island of Santorini in the Aegean Sea (Hübscher et al., 2015;Nomikou et al., 2016b). The wider geological setting includes the Anydros Basin -a 320 region known to be shaped by extensive volcanism resulting in pronounced structural complexity. Owing to past and ongoing tectonic processes, the upper crust is disrupted by major fault and fracture networks. It is dominated by the Kolumbo submarine  volcano, whose activity might have triggered devastating Tsunamis in the past (Nomikou et al., 2016a). Despite the fact that only a single channel was available, the data can be considered of reasonably high-quality. Owing to a short shot interval, the near-offset data set provides dense spatial sampling, which is deemed ideal for high-resolution diffraction imaging. Figure   325 8 shows the single-channel data before processing together with the reconstruction based on coherent diffraction imaging.
Captured is a larger sedimentary basin near the flank of the Kolumbo submarine volcano.
In its unprocessed form the dataset is dominated by reflected energy and only occasionally, small-scale structural complexity is indicated by faint interference patterns. In strong analogy to the marine seismic multi-channel dataset discussed before, the diffraction image appears contrast-rich and highly resolved predominantly in regions where dynamic processes were at work,

330
whereas sedimentary reflections are fully transparent in the reconstruction. Similar to the previous example the dominant features are connected to a major fault system in the right part of the sedimentary basin. This intricate network is thought to represent a major flower structure, which is connected to the rifting regime it is embedded in (Hübscher et al., 2015). Likewise well-resolved are erosional unconformities and the sediment-crystalline-basement interface which are of rugged character with small-scale lateral complexity. To closer inspect the success of the coherent diffraction imaging workflow, Figure 9 compares a 335 close-up (indicated by the frame in Figure 8) of the reflection-dominated input data, the result of diffraction separation and the diffraction-based reconstruction. Quasi-parallel faults are individually recovered with a lateral separation as small as 200 m, which is broadly in the range of the predominant seismic wavelength. To further illustrate the lateral resolution achievable with diffraction imaging, Figure 10 compares an excerpt of a Kirchhoff-migrated image with the corresponding diffraction-based reconstruction. In the latter, aside from the contrast-rich detection of an internal unconformity and the top of the crystalline 340 basement, small-scale sub-parallel faults are recovered at lateral distances above 4 km. Again, the achieved lateral resolution approaches the order of the seismic wavelength, thereby highlighting the high-resolution potential of diffraction imaging. For a more detailed geological interpretation and an example of diffraction-based velocity model building in depth, we refer the reader to Preine et al. (2020).

345
Complementing the two marine seismic examples, here we will briefly demonstrate the successful application of 3D diffraction imaging on land. In 2014 Bob Hardage and Scott Tinker of the University of Texas at Austin decided to make the Stratton 3D dataset, consisting of a migrated reference volume, unprocessed prestack gathers, vertical seismic profiling data, well log information, and other related resources freely available to the research community. The multi-channel seismic data were recorded in four overlapping swaths, which each consisted of six receiver lines separated by approximately 400 m. The acquisition 350 was intended to capture a portion of the Stratton gas field located in Southern Texas, with the aim of better understanding the internal architecture of complicated oil and gas reservoirs. As the gas field is of fluvial origin, dominant structural features include overlapping fluvial channels (Hardage et al., 1994). To illustrate the potential added value of coherent diffraction imaging for 3D seismic interpretation, Figure 11 shows an exemplary migration slice at reservoir level and compares it with an augmented image consisting of the same migration slice overlain by a semitransparent version of the coherent diffraction 355 image. As expected from a thin-bed fluvial reservoir system, the diffraction-based reconstruction indicates a high degree of structural complexity. A comparison of the migration with the augmented diffraction image reveals several spatially coherent diffractive corridors that follow established channel trends and might also be indicative of minor faulting, which is of potential significance for the monitoring of fluid flow.

360
The fourth and final example consists of a ground-penetrating radar (GPR) line that was acquired by the US Geological Survey in the frame of a multi-disciplinary effort to study the impact of seasonal tropical storms on coastal change as part of the Barrier Island and Estuarine Wetland Physical Change Assessment (Zaremba et al., 2016). The campaign was a response to Hurricane Sandy which approached the US East coast in October 2012. A connected storm surge and wave activity caused major alter- ations of the shore line, resulting in a modified coastal topography, geology, and hydrology with immediate impact on regional 365 ecosystems. The GPR measurement campaign was carried out roughly three years after Hurricane Sandy hit the shore and, in conjunction with remote sensing and sedimentological observations, had the main goal of systematically assessing physical changes of the coast induced by seasonal storms in order to update systemic models to improve predictability (Plant et al., 2018). In Figure 12 the exemplary line of the GPR data volume is displayed with the result of diffraction separation and the outcome of coherent diffraction focusing. Along the line, pervasive diffraction can be observed which indicates a strong degree 370 of structural complexity near the surface. At the centre of the line a V-shaped structure is revealed to be responsible for pronounced electromagnetic scattering, suggesting the presence of a strong material contrast. Similar to the single-channel seismic example, owing to the fact that emitting and receiving antennae typically coincide, reflections in GPR data are predominantly sensitive to vertical structural changes, whereas lateral information stemming from channels and other dynamically important erosional structures is largely encoded in the diffracted wavefield.

Discussion and outlook
Aside from arriving at high-resolution structural subsurface reconstructions, diffractions also provide unique illumination in various important data sub-domains, like the zero-offset (single-channel), common-source, or common-receiver configuration.
Illustrated by the marine single-channel example (Figures 8-10), the 3D reconstruction based on reduced land data (Figure 11), and the zero-offset GPR application (Figure 12), diffractors can indeed be viewed as structure-related passive sources, which suggests a systematic use of cost-effective, reduced acquisitions in seismic investigations (Schwarz, 2019a). In addition to the potential of lowered acquisition costs, the strong similarity of diffracted and passive events also suggests continuing technological transfer between controlled-source and passive-source seismology (Li et al., 2020). The introduction of interferencesensitive data abstractions like the considered sign-sensitive n-th root was shown to have a stabilizing effect on diffraction imaging and is also expected to benefit automated coherence analysis as a whole.

Potential and extension of the method
Although not considered here, it should be stated that, similar to the employment of ray tracing or eikonal solvers to arrive at more accurate focusing trajectories in laterally complex media, the projection recipe can likewise be extended to account for more complicated and demanding scenarios. In fact, wavefront tomography builds on the same attribute fields and coherence maps as were fed to the analytic mapping equations (8), but makes use of dynamic ray tracing to perform the subsequent 390 projection into image space. Also, the method was shown to reliably utilize diffracted energy to arrive at a resolved estimate of scatterer locations in depth (Duveneck, 2004;Bauer et al., 2017). While the projection step, like in focusing, could be performed for a pre-defined velocity distribution, projection was shown to lend itself well for a largely automated reconstruction of not only the scatterer locations, but also the macro-velocity structure in an iterative process. As was demonstrated e.g. by Bauer et al. (2019), data-centric mapping lets one tag and track the contribution of each data point into the image, which 395 provides a powerful interface to machine learning techniques, commonly used for image segmentation tasks (Shustak and Landa, 2018). In conjunction with time-reversal imaging (e.g. Mosk et al., 2012), which can likewise be considered a projection technique, coherence arguments and wavefront tomography were recently demonstrated to form a powerful framework for the joint inversion of passive-source and medium properties (Diekmann et al., 2019). Sufficiently dense spatial sampling provided, coherent diffraction imaging, in particular when phrased as a projection problem, is expected to be applicable to passive seismic 400 data, as the problem strongly resembles diffraction imaging in the common-source domain (Schwarz, 2019a). Differential semblance optimization and, more broadly, migration velocity analysis can be viewed as the focusing-based counterpart of wavefront-tomographic inversion and might likewise be investigated in the future.

Limitations and challenges
It needs to be emphasized that, in particular when only reduced datasets like single-channel volumes are acquired, diffractions 405 bear a distinct advantage over reflections, in that their illumination is encoded in various sub-domains of the multi-channel response (Schwarz, 2019a;Preine et al., 2020). However, it must also be noted that the process of diffraction is inherently three-dimensional, which can cause out-of-plane energy to contaminate the data with the potential to result in the occurrence of artefacts in the reconstruction. Provided accurate velocity information is available, out-of-plane contributions in twodimensional (2D) surveys are naturally suppressed, if the scattering structure is located reasonably far off the acquisition plane.

410
However, less distant out-of-plane scattering can hardly be distinguished from valuable in-plane contributions, which is why 2D diffraction-based reconstructions must generally be assessed with care. This is particularly true for single-channel data, for which a reliable velocity model might not be available. In order to gain trust in diffraction-derived velocity information and coherent diffraction images the mere quality of focusing might be complemented by a joint assessment of the reflected wavefield.
Powerful and reliable reflection-based velocity inversion schemes exist and can be used, if sufficient offset information is avail-415 able. Thus, because reflected energy is less likely to stem from out-of-plane structures, the integrated interpretation of reflection and diffraction images can help to improve velocity models and identify off-plane scattering in 2D surveys (Preine et al., 2020).
All these complications become superfluous for sufficiently dense three-dimensional acquisition strategies, which, therefore, are deemed ideally suited for reliable subsurface imaging with the diffracted wavefield. Although coherence arguments were demonstrated to help with handling data imperfections, the highest-quality reconstructions are expected for sufficiently dense 420 spatial data sampling. With the rise of large-N arrays and, in particular, the emerging new data resource of distributed fibreoptic strain sensing (e.g. Jousset et al., 2018), wavefields start to be acquired quasi-continuously, which is expected to extend the reach of coherence analysis and diffraction imaging.

Geological interpretation
In line with previous works on diffraction imaging (Lowney et al., 2020;Preine et al., 2020), the systematic synthetic inves-425 tigations and, in particular, the range of considered challenging field data examples suggest that diffraction imaging bears the potential to shed a unique light on intricate fault and fracture networks and other dynamically relevant features. As indicated in the generalized workflow shown in Figure 2, coherent diffraction imaging lets us arrive at resolved images of lateral discontinuities, which, in some sense, encode the dynamic history and past and present stress states of the crust. In conjunction with conventional reflection-dominated images resulting from migration, diffraction images were shown to be of immediate 430 and complementary use in interpretation (Burnett et al., 2015;Preine et al., 2020).
To decipher the signatures of faults and fracture systems is of immediate importance for the construction of geodynamic models or the simulation and assessment of fluid flow. The suggested incorporation of coherence arguments in constructing images is expected to aid in this interpretation task and interface well with automation strategies that build on machine learning techniques that have their origin in computer vision (e.g. Wu et al., 2019). Generally it can be argued that the positive-valued 435 character of the coherent reconstructions favours the subsequent application of useful tools from image processing. The nonnormalized beam energy (n=1) directly relates to the diffraction's focusing intensity, which is proportional to the square of the beam's amplitude and, therefore, to the strength of the impedance contrast at which diffraction occurred. On the contrary, higher-order versions of the beam energy (n>1) no longer deal with accurate, but rather distorted amplitude and phase information and, accordingly, cannot be used for quantitative interpretation. The same holds for the semblance norm in general, as 440 its intrinsic normalization evens out amplitude discrepancies due to material contrasts of different strength. While all of the coherence measures suffer from the loss of phase information in the final reconstruction, the semblance coefficient, due to its normalization, can be used as a reliable weight for artefact and noise suppression in conventional wavefield focusing. The resulting images have higher quality yet largely preserve amplitude and phase information critical for quantitative geological interpretation.

445
While the presented workflow discusses the best use of the physical information content of the recorded data through diffraction-targeted processing, structural interpretation makes additional use of the growing amount of seismic attributes (Chopra and Marfurt, 2005;Barnes, 2016) -an integral approach of seismic interpretation aiming at mapping geological features. Like coherence (gained via cross-correlation of neighbouring traces in the reflection-dominated migrated image), these attributes are often used on their own to improve the interpretation of fault structures. Alternatively, attributes can be 450 assessed in combination or help in establishing cross-plotting maps (e.g. Endres et al., 2008;Lohr et al., 2008;Torrado et al., 2014;Wang et al., 2015). In this frame, coherent diffraction images can be viewed as physics-guided feature maps that naturally complement more conventional attributes commonly used for interpretation. To additionally foster the bridging from faults to fractures, data acquisition can likewise play an important role (see concept and example in Krawczyk et al., 2015). In nearsurface applications in the field, using shear waves instead of compressional waves for seismic surveying has proven a powerful 455 strategy for increasing structural resolution (e.g. Krawczyk et al., 2012;Beilecke et al., 2016). A combination of the proposed high-resolution imaging workflow with these new forms of data acquisition is expected to shed additional light on subsurface pathways, fault extent and fault connections in the subsurface, which are increasingly important for the assessment of structural integrity and fault behaviour or, ultimately, deformation monitoring in an area.

460
We have presented a simple yet powerful framework to arrive at highly resolved structural images of the upper crust by making use of the diffracted component of the wavefield. By means of controlled synthetic test cases, we introduced and systematically investigated four positive-valued coherence measures which find strong correspondence in visual perception. Based on the prerequisite of reasonably dense spatial sampling, we suggested a generalized workflow for diffraction imaging, in which image formation is either achieved by focusing or by projection of coherent contributions. While synthetic tests suggest the overall 465 robustness of coherence-based diffraction imaging, the investigation of seismic and ground-penetrating radar data acquired on land and in the marine environment emphasize the applicability and complementary nature of diffraction imaging for a broad range of geophysical applications. Owing to its high-resolution potential, the presented workflow helped to delineate smallscale structural features such as fluvial channels, erosional unconformities, and intricate fault and fracture systems, which remain challenging to image by conventional means.

470
Code and data availability. The two suggested strategies for coherent diffraction imaging were implemented in the emergent high-level language Julia. Both, the codes and intermediate and final results are available upon reasonable request.
Author contributions. BS performed the computations and is responsible for the theoretical aspects of the paper. CMK helped strengthening the overall scope and added to interpretational aspects and the discussion of the presented results. Both authors contributed equally to proofreading and the final preparation of the manuscript. the public. Wavefield simulations for the synthetic examples as well as pre-processing and visualisations were carried out with Seismic Unix, whose contributors we wish to acknowledge. We also warmly thank Alexander Bauer and Jonas Preine for continued stimulating discussions.
The thoughtful comments of two anonymous reviewers helped to improve the manuscript.