P Wave Anisotropy Caused by Partial Eclogitization of Descending Crust Demonstrated by Modeling Effective Petrophysical Properties

Seismological studies of large‐scale processes at convergent plate boundaries typically probe lower crustal structures with wavelengths of several kilometers, whereas field‐based studies typically sample the resulting structures at a much smaller scale. To bridge this gap between scales, we derive effective petrophysical properties on the 20‐m, 100‐m, and kilometer scales based on numerical modeling with the finite element method. Geometries representative of eclogitization of crustal material are extracted from the partially eclogitized exposures on Holsnøy (Norway). We find that the P wave velocity is controlled by the properties of the lithologies rather than their geometric arrangement. P wave anisotropy, however, is dependent on the fabric orientation of the associated rocks, as fabric variations cause changes in the orientation of the initial anisotropy. As a result, different structural associations can result in effective anisotropies ranging from ~0–4% for eclogites not associated with ductile deformation to up to 8% for those formed during ductile deformation. For the kilometer‐scale structures, a scale that in principle can be resolved by seismological studies, we obtained P wave velocities between 7.7 and 8.0 km s−1. The effective P wave anisotropy on the kilometer scale is ~3–4% and thus may explain the backazimuthal dependence of seismological images of, for example, the Indian lower crust currently underthrusting beneath the Himalaya. These results imply that seismic anisotropy could be the key to visualize structures in active subduction and collision zones that are currently invisible to geophysical methods and thus can be used to unravel the underlying processes active at depth.


Introduction
Convergent plate boundaries are among the most important sites of crustal reorganization and element recycling. There, crustal material is buried to great depths, recycled into the mantle, integrated into orogenic roots, and in some cases also exhumed back to the surface. All of these processes result in the modification of crustal rocks through metamorphism and brittle and/or ductile deformation. However, this occurs at depths inaccessible to direct observation. Thus, such structures are either studied by geophysical imaging methods or by investigating exhumed rocks that have been metamorphosed and/or deformed in the past (e.g., Austrheim, 1987;Rondenay et al., 2008). Field-based studies of deep processes are restricted to rare exposures where mineral assemblages and structures are not substantially overprinted during exhumation (e.g., Austrheim, 1987;John & Schenk, 2003). In order to properly interpret seismic velocities and deduce the ongoing metamorphic processes associated with large-scale tectonics, we require knowledge of how seismic properties change with depth and lithology (e.g., Kind et al., 2012;Rondenay et al., 2008).
While field-based studies include information down to the micron scale, geophysical imaging techniques employ wavelengths that are only sensitive to kilometer-scale structures (e.g., Bloch et al., 2018;Kim et al., 2019). In addition, the resolution of geophysical imaging is often further limited by the available station coverage and distribution of signal sources. This creates a large gap between the scale at which we image structures with geophysical methods and the scale at which we can observe structures in the field. Subsequently, seismic velocities that are measured in the laboratory or calculated for individual samples may not be representative of the properties of lithological and structural associations on a larger scale. As these structures are smaller than the resolution of seismological methods, the properties of the different constituents will act together as one effective medium (e.g., Backus, 1962;Hudson, 1981;Okaya et al., 2019).
Specifically, eclogitization processes occurring at depth remain difficult to assess, although they are suspected to play a major role in geodynamic processes (Austrheim, 1991;Dewey et al., 1993;Yamato et al., 2019). Eclogitization causes a density increase of crustal material that decreases buoyancy forces and significantly adds to driving forces (e.g., slab pull) at convergent plate boundaries (e.g., Hetényi et al., 2007;Klemd et al., 2011). However, the same density increase also significantly complicates the detection of eclogites at depth as it is combined with an increase of the elastic moduli of the rock. Subsequently, the resulting seismic properties of eclogites become similar to those of mantle peridotites. This makes distinction between the mantle and crust at depth difficult (e.g., Bostock, 2013;Hetényi et al., 2007;Rondenay et al., 2008;Yuan et al., 2000). Nevertheless, partially eclogitized material within a subducting slab shows a range of geometric configurations and orientations of anisotropy in the constituent lithologies, depending on conditions during formation (John & Schenk, 2003;Scambelluri et al., 1995;. It is therefore not necessarily straightforward to transform a measured velocity into a degree of eclogitization.
Eclogites formed from dismembered parts of the subducting crust, for example, at the plate interface, often occur as undeformed boudins in a weaker matrix, typically composed of metasediments (e.g., Hetzel et al., 1998;Pleuger et al., 2005) or serpentinites (e.g., Scambelluri et al., 1995). On the other hand, field-based studies have shown that intraslab eclogitization of crustal rocks is often associated with fluid availability that enhances mineral reactions and ductile deformation, first forming centimeter-thick shear zones (Austrheim, 1987;John & Schenk, 2003). As eclogitization and deformation progress, such shear zones can widen and connect into larger shear zone networks surrounding low-strain domains (e.g., Jolivet et al., 2005) and the shear zones can ultimately reach a thickness of a few hundred meters (Angiboust et al., 2011;Boundy et al., 1997;Raimbourg et al., 2005;). In exposed examples of coherent pieces of partially eclogitized crust, the preserved shear zones rarely reach scales that can be resolved with geophysical methods and the complex associations would thus act as an effective medium at depth (e.g., . In contrast, geophysical imaging methods are used to study large-scale processes active at great depth in collision and subduction zones (e.g., Halpaap et al., 2018). To unravel structures caused by metamorphism coeval with deformation, the receiver function method is of specific interest. It is based on the conversion of P to S waves and vice versa at boundaries with contrasting impedance and therefore mostly sensitive to structural boundaries (Kind et al., 2012). For example, Schneider et al. (2013) imaged a low-velocity zone below the Pamir corresponding to the subducting lower continental crust of the Eurasian Plate. The velocity contrast of this zone with respect to the surrounding mantle, however, decreases below a depth of 100 km, suggesting eclogitization of the downgoing crust. Nabelek et al. (2009) and Schulte-Pelkum et al. (2005) observed a backazimuthal dependence of the retrieved signal in the lower crust of India beneath the Himalaya that suggests a significant large-scale anisotropic fabric within the lower continental crust of India.
Direct estimates of seismic velocities are usually derived from samples that are only a few centimeters in size (e.g., Kern et al., 1996), and extrapolation to scales that are resolvable using geophysical methods relies on poorly supported assumptions, mainly that the composition of the samples is representative of the crust at geophysically relevant scales and that the large-scale organization of lithologies has no relevance. Voigt-Reuss-Hill (VRH) averaging is the standard method to calculate velocities within a medium based on the abundance of individual mineral phases resulting in an average (isotropic) seismic velocity (Hill, 1952). The classic Backus averaging allows calculation of the effective anisotropy of a finely layered medium. It is valid under the assumption that the thickness of individual layers is far smaller than the seismic wavelength (Backus, 1962). Although such averaging schemes are widely used to constrain seismic velocities of various rocks, their capabilities are limited because they are only valid for simple geometries that generally do not capture the structural complexity of real rocks.
To assess these simplifying assumptions, it is necessary to utilize a more sophisticated approach. As a first step in this direction, we focus on the calculation of effective P wave velocities of eclogite-facies associations using a technique based on stress calculations, for a variety of representative geometries. The simplified geometries are derived from field observations on the island of Holsnøy in the Bergen Arcs (Norway), where a >70 km 2 large complex of partially eclogitized lower continental crust that provides an excellent coherent laboratory to study the geometries that are established during eclogitization is exposed.

Geological Setting
The exposed lower continental crust on the island of Holsnøy (Bergen Arcs, western Norway) has been partially eclogitized during the Caledonian orogeny (Austrheim, 1991). The rocks belong to the Lindås nappe, which, together with the Dalsfjord and Jotun nappe complexes, represents the lower crust of the former Jotun microcontinent that constituted part of the pre-Caledonian hyperextended margin of Baltica (Andersen et al., 2012;Jakob et al., 2019). The Lindås nappe is for a large part composed of anorthositic granulites that experienced Proterozoic granulite-facies P-T conditions of~1 GPa and~800°C, at~950 Ma (Austrheim & Griffin, 1985). The P-T conditions in the following~500 Myr are unclear. The rocks, however, show no signs of significant alteration before the Scandian Caledonian collision and likely cooled to conditions reflecting mid to lower crustal conditions (Jamtveit et al., 1990).
During the Caledonian collision, the Jotun microcontinent constituted the leading edge of Baltica, which was integrated into the collision wedge as the lower plate (Corfu et al., 2014). Subsequently, the Lindås nappe was subjected to peak eclogite-facies conditions of~2 GPa and~750°C at 429 Ma (Bhowany et al., 2018;Glodny et al., 2008;Jamtveit et al., 1990;Zhong et al., 2019). Large volumes of the dry granulite-facies rocks, however, remained metastable and were thus preserved (Austrheim, 1987;Jackson et al., 2004). Eclogitization is linked to fluid availability and was facilitated along shear zones but also progressed into the rock volume as a static overprint (Austrheim, 1987;. Fluid infiltration was likely initiated via brittle fractures, which provided fluid pathways within an otherwise dry rock (Austrheim, 1990;Jamtveit et al., 1990).
This heterogeneously distributed transformation resulted in a complex mixture of eclogites and granulites ( Figure 1). The resulting lithologies can be divided into six categories based on the abundance of eclogite and the associated structural relationships (Boundy et al., 1992;. Next to the mostly unaltered granulite (<20% eclogite), small-scale eclogitization features are distinguished into granulites cut by eclogite-facies shear zones a few centimeters wide and granulites with eclogitized patches that are not associated with ductile deformation (both 20-50% eclogite). With progressive eclogitization, these evolve into the so-called eclogite breccia, which can be described by two endmembers: sheared eclogite breccia composed of a strongly sheared eclogite matrix containing preserved granulite blocks and unsheared eclogite breccia, where the eclogite matrix was not subjected to pervasive ductile deformation (50-90% eclogite). Ultimately, shear zones that are up to a few hundred meters thick and are almost entirely composed of eclogite with little to no preserved granulite evolve (>90% eclogite).

Finite Element Calculations
The aim of this study is to obtain effective P wave velocities and the corresponding P wave anisotropy from variably eclogitized lower crustal rocks based on observed 2D geometric arrangements that act as an effective medium. The information extracted in the field is simplified and translated into numerical models with the goal of capturing the most essential properties of the observed field relationships. As the field observations along approximately planar exposures would require the introduction of essentially arbitrary additional parameters to create plausible 3D models, this contribution focuses on 2D numerical modeling. Whereas this approach inevitably results in some differences in the inferred seismic velocities, compared with what would be obtained for the unknown true 3D structure, the 2D models nevertheless are expected to provide a good approximation and therefore enhance our understanding of how geometries affect petrophysical properties from the outcrop to the map scale. It has to be noted that here we focus on P waves, as in our formulation, the out-of-plane polarization of S waves is not included.
Both the effective medium and the individual rock types are treated as linear elastic anisotropic material for which Hooke's law gives the relationship between stress (σ ij ) and strain (ε kl ): σ ij ¼c ijkl ε kl : In 2D, the 2 × 2 × 2 × 2 elastic tensor, which we represent by a symmetric 3 × 3 matrix in Voigt notation (using the mapping 11 → 1, 22 → 2, and 12 → 3), is sufficient to fully describe the in-plane anisotropy: Due to symmetry considerations, c 13 , c 23 , c 31 , and c 32 are expected to be zero, and c 12 should be equal to c 21 . One way of obtaining the effective properties is to run numerical experiments solving the elastodynamic wave equations and recording the time necessary for a wave to travel through the medium (e.g., Saenger et al., 2004). Other studies have used asymptotic expansion homogenization to calculate the effect of mineral orientation on seismic anisotropy solving for the 6 × 6 elastic tensor in 3D (Naus-Thijssen, Vel et al., 2016). Alternatively, we calculate the P wave velocities from the elastic tensor of the effective medium using the formulas for transversely isotropic media (Mavko et al., 2009): The equations fully describe the P wave velocity in 2D because our calculations assume no out-of-plane properties, that is, plane strain. Thus, the reduction of a transversely isotropic medium to 2D, with the symmetry axis within the plane, results in an anisotropic medium and is sufficient to describe the anisotropic elastic properties in 2D. The individual components of the 2D elastic tensor (c ijkl ) of the effective medium are calculated from the stresses and strains calculated in a set of numerical experiments. For this purpose, three experiments ( Figure 2) are performed for each geometric model, applying different boundary conditions: (1) the area of interest is compressed along the y axis along the upper and lower boundaries by imposing a fixed displacement. Along the left and right boundaries, displacement in the x direction is zero.
(2) The medium is compressed horizontally, that is, along the x axis. In this case, displacement in the y direction is zero along the top and bottom boundaries.
(3) Finally, simple shear is enforced along the top and bottom boundaries, that is, displacement to the right along the top boundary and to the left at the bottom boundary, resulting in shear parallel to the x axis. A fourth experiment (simple shear parallel to the y axis) was used for validation and yielded the same results as Experiment 3, as is required from the symmetry of the elasticity tensor.
The three experiments result in a set of nine equations for six unknown components of the stress tensor, so only six of these equations are needed. Due to the setup of each experiment, specific strains are zero, which allows the simplification of the equations to for Experiment 1; for Experiment 2; and for Experiment 3.
To extract the elastic properties of the effective medium, strain (ε kl ) and stress (σ ij ) are averaged across the domain. The boundaries at which the displacement for each of the experiments is enforced are kept far away from the medium of interest to avoid boundary effects. Strain and stress are then averaged only across the domain which constitutes the medium of interest (red in Figure 2). In the 20-m scale models (Figures 3a and 3b), both lithologies are in contact with the inner boundary in some cases.
Here the eclogite was extended into the area between the inner and outer boundaries to avoid edge effects.  (a) Small-scale eclogite-facies shear zones representative of an area of~20 × 20 m. The example shown here contains 30% eclogite. For the calculations with other eclogite abundances, the thickness of the shear zone was varied accordingly. (b) Small-scale static eclogite overprint representative of an area of~20 × 20 m. The example shown here contains~30% eclogite. For calculations with other eclogite abundances, the size of the eclogite patches was varied accordingly. (c) Sheared eclogite breccia with regularly oriented granulite blocks. The example is representative of an area of~100 × 100 m and~70% eclogite. The size of the granulite blocks remains the same throughout all calculations. To perform calculations with different eclogite abundances, the abundance of granulite blocks was altered. (d) Unsheared eclogite breccia with the same variations as in (c). Below each image, the corresponding properties of eclogite and granulite used for the calculations are given. Each column represents one model series. The percentage gives the strength of the P wave anisotropy of the corresponding rock, and the arrow gives the orientation of the fast P wave direction used for the calculations. L and H indicate whether the higher-or lower-anisotropy version was used.
The P wave velocities of the effective medium can then be calculated from the resulting 2D elastic tensor. Bulk density is obtained by calculating the mean of the densities weighted by the area of the granulite and eclogite used for the calculation.
The calculations are performed using the Galerkin finite element method (FEM) employing an irregular triangular grid. Meshing is done with the mesh generator triangle (Shewchuk, 1996). Each triangular element consists of six nodes in which the displacement field is calculated with quadratic interpolation and Gauss quadrature integration at three points.
The method of obtaining the P wave velocity described above was tested and benchmarked using an anisotropic layered medium, a problem for which an analytical solution exists (Backus, 1962). Benchmarking was performed on a regular grid and reproduced the analytical solution within machine level precision.

Properties of the Implemented Lithologies
The physical properties for each element representing the different materials are given by the elastic tensor of the corresponding lithology, that is, granulite or eclogite. Representative elastic tensors were calculated from the velocity measurements (x-z plane) in . Those measurements were performed on a true triaxial multianvil press using the ultrasonic pulse transmission technique (Kern, 1978) with varying pressure and temperature between ambient conditions at 600 MPa and 600°C, respectively. From this data, each component can be calculated separately except for c 12 , as this would require information on the variation of elastic wave speeds along oblique directions not available from laboratory measurements. Therefore, we used the mean P wave velocity between the x and z axes to approximate the velocity of a P wave traveling at a 45°angle to the foliation. This results in an almost elliptical anisotropy with the ellipticity parameter η κ (Brownlee et al., 2017;Kawakatsu, 2016) varying between 0.97 and 1.03. Brownlee et al. (2017) have shown that off-axis anisotropy deviates systematically from elliptical symmetry for rocks with high anisotropy. However, this effect is most pronounced for rocks with high mica and/or quartz contents, which is not the case for the Holsnøy samples shown in . The assumption of near-elliptical anisotropy made here is thus in agreement with the scaling laws of the off-axis anisotropy proposed by Brownlee et al. (2017). In order to test the relative influence of the intrinsic properties of the constituting lithologies and the geometries themselves, we used two different eclogites and two different granulites (Table 1). Because the eclogites measured by  were all collected from the main shear zones exposed on Holsnøy, they all have a high P wave anisotropy. In order to estimate effective properties for statically eclogitized areas, where the eclogite would likely have a lower initial anisotropy, we assumed a lower velocity in the x direction for one of the samples (N-101 in , thus giving a lower P wave anisotropy of 4%, which is in accordance with others reported from Holsnøy (Fountain et al., 1994). Specifically, we chose to use the velocity measured at lower confining pressure (600 MPa). This way, while the velocity is artificially reduced, it is still a function of the existing mineral assemblage.
The calculations feature four different categories of structural associations as they are evident from the field (e.g., Austrheim, 1987;Raimbourg et al., 2005 Figure 3d). For each of the categories, a series of calculations was performed, systematically varying the main configurations that can be observed in the field. These are (a) abundance of eclogite (10-50% eclogite for 20-m scale and 50-90% eclogite for 100-m scale), (b) orientation of the main foliation of the constituting lithologies (Figures 3 and 4), and (c) strength of the deformation fabric in the lithologies (static vs. dynamic eclogitization; . The orientation of the foliation used for the calculations is given in Figure 3 (as sketches) and Figure 4, where XX means that the fast axis and thus the foliation of the two lithologies are parallel and XY means that they are perpendicular to each other. The strength of the deformation fabric is also given in Figure 3 (in % anisotropy) and Figure 4 with the following notation: E,H-higher-anisotropy eclogite; E,L-lower-anisotropy eclogite; G,H-higher-anisotropy granulite; and G,L-lower anisotropy granulite. Note. The velocities (V P and V S ), densities, and anisotropies were taken from . The star indicates that V PX of N-101 was adjusted so that an anisotropy of 4% results (see text). Anisotropy was calculated as 100*(V PX -V PY )/V Pmean . Velocities (V) are given in km s −1 , density (ρ) in kg m −3 , and anisotropy (A) in %.
The choice of orientation of these lithologies is related to their structural setting. For the small-scale eclogite shear zones, the fast axis of the eclogite is oriented parallel to the shear zone and its foliation as would be established during ductile shear (e.g., Bascou et al., 2001). The fast axis of the granulite is oriented either parallel or perpendicular to the shear zone, thus representing a situation where the shear zones are established parallel to the foliation of the granulite or perpendicular to it. On Holsnøy, both of these scenarios are present in the field; however, most shear zones develop obliquely to the granulite-facies foliation. The two scenarios introduced into the calculations are thus endmember representations of the chosen field example.
Patches of statically eclogitized material on the other hand typically develop parallel to the granulite foliation. For those calculations, the fast axis of the granulite and subsequently the foliation are horizontal. The fast axis of the eclogite is either parallel to the granulite foliation as typically observed in the field or perpendicular to it. We chose to include this case as field observations suggest that static eclogitization features can also crosscut the granulite foliation, specifically in cases where the abundance of eclogite is large enough for individual granulite blocks to start moving independently. This setup thus also represents endmember scenarios.  (b) is to the right of (b) and the legend for (c) and (d) is to the right of (d). Each model series is shown by two connecting dashed lines. The upper line represents the maximum velocities, and the lower line represents the minimum velocities. The legend is given in the following scheme: XX-foliation and fast axis of granulite and eclogite are parallel; XY-foliation and fast axis of granulite and eclogite are perpendicular (shown as sketches in Figure 3); E,H or E,L-eclogite with high anisotropy or low anisotropy, respectively; and G,H or G,L-granulite with high or low anisotropy, respectively. Higher and lower anisotropy is indicated in % in Figure 3. Solid lines indicate Voigt-Reuss-Hill averages in the same color scheme as the modeling results averaged from the elastic tensor implemented in the corresponding calculations. Shown are calculated velocities parallel to x and y.
Finally, for the 100-m scale calculations of sheared and unsheared eclogite breccia, the eclogite foliation is horizontal in the models. For the sheared eclogite breccia, this is straightforward, representing the foliation formed during ductile shear. In the case of unsheared eclogite breccia, this was done to make comparison between the two easier. The granulite blocks in the sheared case are often aligned nearly parallel to the eclogite foliation with their preserved internal (granulite) foliation. However, rotation of these blocks is not always complete, and the granulite foliation can appear random in some cases. For the unsheared eclogite breccia, the granulite blocks in the field can indeed have all possible orientations. Subsequently, the implemented orientations of the fast axis of granulite and eclogite also cover the endmember cases of what can be observed in the field.

P Wave Velocity and Anisotropy of Small-Scale Eclogitization Features
The process of eclogitization, as it can be studied on Holsnøy, is driven by two contrasting endmember mechanisms: eclogitization proceeding along shear zones (Figure 3a) or developing as a static overprint (Figure 3b). In the first type, eclogitization proceeds along shear zones that widen progressively with time, and a wide variety of shear zone thicknesses is found throughout the field area (Austrheim, 1987). Hence, we calculated P wave velocities for 20 examples with varying shear zone thickness as well as varying elastic properties of the eclogite and granulite implemented in the models (Figure 3, Table 2).
Comparing all models shows that the calculated P wave velocities with higher-anisotropy (stronger deformation fabric) granulites are, in general, higher than those with lower-anisotropy granulites. Furthermore, with increasing shear zone thickness (i.e., amount of eclogite), the P wave velocity in both the slow and the fast P wave directions increases linearly (Table 2; Figure 4a). The one exception to this trend is given by the models that feature both higher-anisotropy eclogite and the higher-anisotropy granulite, with the fast axis of both rocks oriented perpendicular to each other: the fast axis of the eclogite parallel to the shear zones and the Note. Velocities are given in km s −1 and anisotropy is given in %. For each model, the properties of the granulite and eclogite used for the calculation are indicated with the following scheme; L: low anisotropy, H: high anisotropy, X: fast axis is oriented horizontally, and Y: fast axis is oriented vertically.
fast axis of the granulite perpendicular to them. For this geometry, the resulting P wave velocities for the fast and slow axis of the effective medium converge up to an eclogite abundance of~30% and then diverge toward higher eclogite abundance. This coincides with a change of the orientations of the fast and slow directions. In this scenario, the velocity perpendicular to the shear zones is almost constant. The velocity parallel to the shear zones, however, increases significantly from 7.58 to 7.96 km s −1 with increasing eclogite abundance. The fast axis is thus perpendicular to the shear zones for an eclogite abundance <30% and parallel to the shear zones from~30-50% eclogite abundance. The orientation of the slow direction rotates progressively away from the trend of the shear zones, that is, toward the slow direction of the eclogite. In all other model sequences, the fast direction is parallel to the shear zones and the slow direction is perpendicular.
The corresponding P wave anisotropy also increases with increasing shear zone thickness and reaches a maximum value of 7.1% ( Figure 5). In most models, this increase is near-linear with increasing eclogite abundance. In contrast, the resulting P wave anisotropy of those calculations featuring a higher-anisotropy granulite with the fast axis oriented perpendicular to the shear zone decreases betweeñ 10% and~30% eclogite abundance and then increases until~50% eclogite abundance. Finally, the P wave anisotropy at~50% eclogite abundance returns to approximately the same value of~2-3%, as the P wave anisotropy at~10% eclogite abundance.
In general, the resulting P wave anisotropy is larger when the fast axes of both granulite and eclogite are oriented parallel to the shear zone, compared with those examples where the fast axis of the granulite is

10.1029/2019GC008906
Geochemistry, Geophysics, Geosystems oriented perpendicular to the shear zone. At lower eclogite abundance, the calculations implementing a higher-anisotropy granulite result in a higher anisotropy of the effective medium, while the results at higher eclogite abundance indicate that the anisotropy of the effective medium is higher if the granulite has a lower intrinsic anisotropy ( Figure 5).
Additionally, eclogitization on Holsnøy also proceeded statically without significant ductile deformation (Jamtveit et al., 2000;. Here, eclogitization most commonly advances parallel to the granulite foliation. For this case, we also calculated 20 different examples, varying both the abundance of eclogite and the elastic properties of the granulite and the eclogite (Figure 3b; Table 2).
The resulting P wave velocities show a similar trend as those from the examples featuring small-scale shear zones. Both the velocities of the fast and the slow axes increase linearly with increasing abundance of eclogite ( Figure 4c). Further, the models featuring granulites with higher anisotropy result in faster P wave velocities of the effective medium than the models implementing lower-anisotropy granulites. Additionally, the P wave velocities are in the same range as the ones calculated for small-scale shear zones.
The orientations of the fast and slow axes are typically constant with the fast axis being parallel to the granulite foliation (horizontal) and the slow axis perpendicular. Only in the calculations where the lower-anisotropy eclogite and granulite with the orientation of the fast axes perpendicular to each other are implemented does the resulting orientation change slightly. These calculations indicate that the fast axis remains horizontal (i.e., parallel to the granulite foliation) while the slow axis rotates slightly away from the initial vertical orientation.
The P wave anisotropy shows a variable trend comparing the different models. The medium featuring the higher-anisotropy granulite, with both the fast axes of the granulite and eclogite oriented parallel to each other, results in a P wave anisotropy of~4% that essentially does not change with varying eclogite abundance (Figure 5d). The same is observed for the models featuring the lower-anisotropy granulite with the fast axes of the granulite and eclogite being oriented perpendicular to each other. Here, the resulting P wave velocity remains relatively constant around 1-2%.
In contrast, the resulting anisotropy of the other two model types changes with increasing eclogite abundance. The sequence featuring a lower-anisotropy granulite with the fast axis oriented parallel to the fast axis of the eclogite increases from~2% to~4%, while the sequence featuring the higher-anisotropy granulite with the fast axes of the two rocks oriented perpendicular to each other decreases from~3% to <1%.

P Wave Velocity and Anisotropy of the Eclogite Breccia
With increasing degree of eclogitization, the so-called eclogite breccia develops, which is composed of an eclogite matrix that surrounds preserved blocks of granulite (Boundy et al., 1992). On Holsnøy, the eclogite breccia can be divided into two endmember types : the sheared eclogite breccia is characterized by a strongly sheared and foliated eclogite matrix, while the matrix of the unsheared eclogite breccia is diffuse and less strongly foliated.
We calculated 20 examples for each of the two types, varying the abundance of eclogite and the elastic properties of the granulite and eclogite (Figures 3c and 3d; Table 2). For all examples of the sheared eclogite breccia, the P wave velocities increase linearly with increasing eclogite abundance. All fast axes and all slow axes converge toward higher eclogite abundances, thus giving fairly distinct maximum and minimum P wave velocities at high eclogite abundances that are independent of the elastic properties of the granulite implemented in the model. (Figures 4 and 5). The slope of the linear increase for the different models is similar to the models dealing with small-scale shear zones. Further, the fast axis of the effective medium in all models is parallel to the shear plane (horizontal) and the slow axis is perpendicular. Additionally, the P wave velocities at~50% eclogite abundance agree well between the models for small-scale shear zones and the sheared eclogite breccia at the same eclogite fraction.
As in the case of the small-scale shear zones, the P wave anisotropy calculated for the sheared eclogite breccia increases nearly linearly with increasing eclogite abundance reaching 7-9% at~90%. Further, P wave anisotropy is consistently higher for models where the fast axes of the granulite and the eclogite are oriented parallel. In this scenario, the anisotropy reaches its maximum when both the granulite and the eclogite have 10.1029/2019GC008906 Geochemistry, Geophysics, Geosystems a high anisotropy. If the fast axis of the granulite, however, is perpendicular to the fast axis of the eclogite, the resulting anisotropy is higher when the implemented granulite has a lower anisotropy.
The P wave velocities calculated for the unsheared eclogite breccia show the same general trends as those for the sheared eclogite breccia (Figure 4d). The only deviation results from the examples implementing granulite and eclogite with their anisotropy perpendicular to each other. Here the calculations result in a change of the orientation at high eclogite abundances.
The trends of the P wave anisotropy of the unsheared eclogite breccia in all calculated examples are lower than the comparable examples of the sheared eclogite breccia ( Figure 5). Most sequences, however, also slightly increase with increasing eclogite abundance, except for those where a lower-anisotropy eclogite is paired with the higher-anisotropy granulite, both of which have their fast axes parallel to each other. In that case, the P wave anisotropy is nearly constant at~3.7% (Table 2).

Discussion
Many studies have calculated or measured P wave velocities of various metamorphic rocks with the aim of interpreting the results of large-scale geophysical imaging techniques (e.g., Almqvist & Mainprice, 2017). However, the sample sizes used for these interpretations are typically far below the resolution of geophysical studies. It is thus essential to understand how geometries formed at depth during ongoing eclogitization shape the seismic properties of the effective medium in combination with the (anisotropic) seismic properties of the constituent rocks.
We distinguish two geometrical contributing factors in order to characterize their influence separately.
(1) The first is the configuration that the different lithologies have to one another on the outcrop scale or larger. This includes, for example, eclogite shear zones that crosscut granulites. In the following, this will be referred to as external geometry, as it involves the relationship of the lithologies to each other but not specifically the properties of the constituting lithologies themselves. (2) The second contributing geometrical factor will be referred to as internal geometry. It highlights the properties of the lithologies themselves by characterizing the relationship between the directional dependence of the elastic properties of the different lithologies that is caused by, for example, crystallographic preferred orientations (CPOs) or shape preferred orientations (SPOs). The internal geometry thus distinguishes whether the fastest velocities of the eclogite and granulite are parallel or oblique to each other.

Effective Properties of 20-and 100-m Scale Structures
Essentially, the P wave velocities calculated for the different geometrical setups show that the velocities are controlled by the velocities of the constituent rocks and their proportions (Figures 4 and 5). This has been accepted and applied by previous studies by calculating, for example, VRH averages (Hill, 1952) and linking those with the CPOs of the mineral phases (e.g., Hacker et al., 2014;Llana-Funez & Brown, 2012;Worthington et al., 2013). Most of these studies, however, obtain information from the thin section scale to recognize crustal-scale processes or to interpret the results from large-scale geophysical imaging studies. The results presented in this study indicate that VRH averages calculated from outcrop-scale features are sufficiently precise to estimate the effective properties on a variety of scales ( Figure 4). Essentially, the external geometries that are representative of eclogitization of crustal rocks have only limited influence on the resulting P wave velocities. Only in isolated cases are the velocities modified, thus deviating from the calculated VRH averages (Figures 4a and 4c). Here a minor geometric effect is plausible; however, this effect results in a maximum modification of <0.2 km s −1 of the P wave velocity and is thus negligible in the context of large-scale crustal processes. In fact, most studies that distinguish between the effect of CPO versus SPO on the thin section scale conclude that the effect of SPO is negligible (e.g., Zhong et al., 2014), although there is evidence that, at least in the deep mantle, SPO does produce seismic anisotropy (Faccenda et al., 2019).
Furthermore, although Holsnøy serves as an example here, this observation can be transferred to other exposures of partial eclogitization and is thus likely representative in more general terms. The shear zones explored here, for example, constitute extreme cases of geometric arrangement as shear zones are well ordered and have a significant lateral extent. It could thus be expected that the influence of, for example, isolated eclogite boudins or blocks as reported from other exposures (e.g., John & Schenk, 2003;Locatelli et al., 2019;Mørk, 1985) on effective P wave velocities would be even smaller.
However, P wave anisotropy varies between the different geometrical configurations ( Figure 5). In this context, our results reveal the importance of the internal geometry compared with that of the external geometry (Figures 4 and 5). As discussed above, the external geometry only has a minor effect on the P wave velocities and anisotropy of the effective medium. The variations of anisotropy for the different configurations tested by us are thus controlled by the internal geometry. The most important factor is the anisotropy of the constituent lithologies that are necessary to produce significant anisotropy of the effective medium. Additionally, the effective anisotropy is strengthened or weakened by the relationship of the individual anisotropies of the lithologies. Anisotropies are higher if the fast axes of the lithologies are aligned but not higher than the highest contributing anisotropy ( Figure 5). Further, our results demonstrate the predominance of the higher-anisotropy lithology. The fast axis of the effective medium is parallel to the anisotropy of the matrix lithology (i.e., in line with the fabric of granulite or eclogite), if the difference in anisotropy between the lithologies is small, or parallel to the higher-anisotropy lithology, even if this lithology is less abundant ( Figure 5). This means that a strongly deformed rock, such as eclogite in shear zones, controls the overall anisotropy even at low abundances. The predominance of the higher-anisotropy phase has also been demonstrated on the rock scale considering, for example, the alignment of mica (e.g., Naus-Thijssen, Goupee, .

Effective Properties on the Kilometer Scale
Combining our results with field observations provides the opportunity to understand how partial eclogitization of crustal rocks alters the seismic properties on a scale significantly larger than what can be measured in the laboratory. Our results suggest that P wave velocities are almost entirely controlled by the velocities and abundances of the constituting rocks ( Figure 4). Essentially, there is no difference in the P wave velocities between rocks that have formed through static eclogitization and those that formed while undergoing ductile deformation. Neither the finite geometries nor the intrinsic seismic anisotropy of the granulites and eclogites have a significant impact on the resulting isotropic average bulk velocities, and the variations that can be distinguished are minor. The P wave anisotropy, however, is influenced strongly by the anisotropy of the rocks that form the effective medium ( Figure 5). Further, our results show that the rock with the higher anisotropy controls bulk anisotropy. In any case, the exemplary geometries discussed above are still far smaller than what can be resolved with large-scale geophysical methods. Therefore, we used these results to extract bulk properties of the effective medium at a scale that could be resolved by large-scale geophysical imaging ( Figure 6). Accordingly, we used an area on Holsnøy that is Geochemistry, Geophysics, Geosystems ~3.9 × 4.6 km in size (Figures 1 and 6a) and provides a coherent natural laboratory for eclogitization-related structures. The geometries are based on the map shown in . As properties for the different map units, we implemented the resulting elastic tensor of the examples shown above, choosing one representative example for each of the geometric configurations, that is, sheared eclogite breccia at~75% eclogite with the fast axis of higher-anisotropy eclogite and higher-anisotropy granulite parallel to each other and unsheared eclogite breccia at~71% eclogite with the fast axis of the lower-anisotropy eclogite and the higher-anisotropy granulite parallel to each other (Table 2). For pure eclogite and granulite, we chose the higher-anisotropy versions ( Table 1) that were also used for the calculations discussed above . The elastic tensors were rotated so that the fast axis is parallel to the structures presented by .
Additionally, we implemented a second (more precise) model that also includes smaller-scale structures (Figure 6b). Here we implemented the small-scale eclogite shear zones at~31% eclogite with the lower-anisotropy granulite and the higher-anisotropy eclogite (Table 2) and the small-scale static eclogite at~32% eclogite with the higher-anisotropy granulite and the lower-anisotropy eclogite ( Table 2).
The resulting P wave velocities of both models are in the range of 8.0 km s −1 (fast axis) to 7.7 km s −1 (slow axis), that is, within the expected range of the measured P wave velocities between granulite and eclogite . Similar velocities are also reported from geophysical studies dealing with active convergent settings, typically in the range of 7-8 km s −1 (e.g., Nabelek et al., 2009;Schulte-Pelkum et al., 2005;Sippl et al., 2013). Additionally, our calculation predicts a P wave anisotropy of 3.9% for the simpler model ( Figure 6a) and 3.3% for the model that includes the small-scale structures (Figure 6b). These values are in the range of what is generally reported from eclogites and granulites (e.g., Brown et al., 2009;Worthington et al., 2013). However, it has to be noted that the anisotropy presented here is representative for the effective medium on a kilometer scale and not only for single (hand specimen-sized) samples.

Implications for Imaging of Continental Collision
The exposures on Holsnøy have been studied extensively because they provide a rare example to study partial eclogitization of a coherent piece of continental crust (e.g., 1987, Bjørnerud et al., 2002Jamtveit et al., 1990Jamtveit et al., , 2000Jolivet et al., 2005). These exposures are also widely considered as representative of how crustal rocks transform to eclogites (e.g., Austrheim, 1990). Specifically, Holsnøy is often regarded as an ideal analog to the underplating crust of India below the Himalaya (e.g., Jamtveit et al., 2019;Labrousse et al., 2010). P wave velocities below the Tibetan Plateau are suggested to be >7.0 km s −1 , which was interpreted to repre-sent~30% eclogitization (Schulte-Pelkum et al., 2005). Our calculation for Holsnøy is representative of~50% eclogitization and yields slightly higher velocities. It thus seems possible to estimate the degree of eclogitization based on P wave velocities. However, this is only feasible if the backazimuthal distribution is sufficiently representative (Nabelek et al., 2009;Schulte-Pelkum et al., 2005).
The retrieved P wave anisotropy of 3-4% from our model is sufficiently high that it could result in a backazimuthal dependence of the retrieved signal in seismological studies. Additionally, our calculations of P wave anisotropy of the different structural associations that could be expected in a partially eclogitized crust show how different geometries can cause high P wave anisotropy ( Figure 5). In the Himalaya-Tibet collision system, where the lower crust of India is imaged below the Himalaya (Jackson et al., 2004;Labrousse et al., 2010), using the receiver function method, it has been shown that the retrieved signal of the Moho is sharp using earthquakes coming from the north, while the Moho cannot be clearly imaged using earthquakes arriving from the south, suggesting an anisotropic fabric within the buried crust (Nabelek et al., 2009;Schulte-Pelkum et al., 2005). Nabelek et al. (2009) propose that this fabric is caused by the imbrication and rotation of a stratified lower crust, excluding eclogites as the cause for the anisotropy because eclogites typically have anisotropies <4%. However, our results show that partial eclogitization of the lower crust does indeed produce considerable anisotropies at the scale sampled by geophysical imaging techniques. Moreover, as shown by our results, the effect of external geometry on seismic anisotropy is limited, suggesting that simple layering or imbrication might not produce sufficient seismic anisotropy on this scale. Our results provide an alternative explanation for the structures observed below the Himalaya. We suggest that considering the P wave velocities reported and the backazimuthal dependence (Nabelek et al., 2009;Schulte-Pelkum et al., 2005), eclogitization of the crust along ductile shear zones, similar to those exposed on Holsnøy, seems the more likely explanation.
Additionally, both kilometer-scale models we present here suggest that the fast axis of the shear zone system is oriented WNW-ESE. At least in a qualitative sense, this suggests that during ongoing eclogitization, when this anisotropy was established, it was dipping toward the upper plate as is also evidenced by the top-east kinematics of the shear zone system Raimbourg et al., 2005). Geophysical imaging suggests a northward dipping fabric within the lower Indian crust (Nabelek et al., 2009;Schulte-Pelkum et al., 2005), that is, dipping toward the Asian plate, consistent with a top to the south shear sense. Our results demonstrate that propagating eclogite-facies shear zones would produce a fabric and subsequent anisotropy with a similar orientation. The scale of those shear zones is actually a minor issue, since our results show the same dependence of effective medium properties on constitutive lithologies independent of the scale.

Implications for Imaging of Oceanic Subduction
Although the rocks on Holsnøy originate from continental crust, some implications for oceanic subduction settings can nevertheless be explored. In many geophysical studies of subducting oceanic plates, the descending crust is clearly imaged at shallow depth but loses its seismic signal at greater depth (e.g., Bostock et al., 2002;Pearce et al., 2012;Rondenay et al., 2008;Yuan et al., 2000). This decrease of the seismic signal is typically interpreted as due to a decreased impedance contrast between descending crust and mantle rocks caused by eclogitization. This is often accompanied by an increase of the dip angle in the Wadati-Benioff zone that indicates a kink in the slab geometry (e.g., Halpaap et al., 2018;Klemd et al., 2011;Yuan et al., 2000). While the subducting crust is invisible to seismological studies at this point, its presence is evidenced by the Wadati-Benioff zone and the inferred kink of the slab has been proposed as a possible geometric obstacle that inhibits exhumation of crustal material subducted beyond that point and is therefore potentially vital to understand subduction zone processes (Klemd et al., 2011). Additionally, kinking on this scale must cause internal deformation of the subducting slab. Whether or not this deformation is localized or homogeneously distributed and how this deformation process affects ongoing eclogitization of the slab is uncertain. However, utilizing seismic anisotropy and the subsequent backazimuthal bias on the retrieved seismic signal might prove to be a powerful tool to unravel these processes in active subduction zones. In this context, although reliable imaging of the crustal anisotropy at these depths is still challenging, seismic anisotropy of the subducted oceanic crust might make it possible to image it to larger depth and illuminate an otherwise invisible slab.

Conclusions
We calculated P wave velocities and the corresponding P wave anisotropy for various geometries, which are representative of partially eclogitized crust. The results show that dynamic eclogitization, associated with shear zone formation, can cause a high P wave anisotropy that increases with increasing eclogitization. The anisotropy of the effective medium is generally controlled by the anisotropy of the matrix or by the contributing lithology that has the highest anisotropy, even if this lithology is less abundant than the other contributors. Consequently, patches of static eclogitization produce a comparatively low P wave anisotropy, which is in some cases independent of the amount of eclogitization. The (external) geometric configuration of the lithologies has little to no effect on the seismic properties of the effective medium.
Our results link partial eclogitization with geophysical observations at active convergent plate boundaries. Previously, significant anisotropy due to eclogitization in deeply buried or subducted crust has been excluded as eclogites are typically not strongly anisotropic. Contrary to this, our results demonstrate that significant anisotropy due to partial eclogitization of crustal material on a kilometer scale is a likely explanation for the discrepancy of the signals retrieved from different backazimuths in seismological studies. For example, the structures seen below the Himalaya are likely anisotropic due to the formation of eclogite-facies shear zones within the lower Indian crust. Additionally, our results strongly encourage the utilization of seismic anisotropy as a tool to visualize the structural associations at depth, thus aiding the extraction of the underlying mechanisms active during ongoing eclogitization of crustal material.