Spacing and Strain During Multiphase Boudinage in 3D

Most models of brittle boudinage predict a dependency of fracture spacing on thickness of the boudinaged layer. This basic relationship can be distorted in the case of multiphase boudinage, where structural inheritance, possibly combined with time evolution of rheology affects boudin geometries, but is not recognized in 2D outcrops. Here we present analyses from a metre-scale, serially sectioned sample of marble from Naxos, Greece, containing a layer of amphibolite, which is known to have experienced two phases of ductile drawn boudinage followed by five generations of brittle boudinage. We show that brittle boudinage depends on the ratio of grain size and layer thickness of the amphibolite. In very thin layers (few grains across), strain is diffuse throughout the entire layer, leading to macroscopically homogeneous stretching. Strain localisation starts when layer thickness exceeds 5 grains: narrow tensile necks and shear zones (35– 45° dip) develop where the layer is <10–20 grains thick. Thicker layers (<30–40 grains) show extensional (mode 1) fractures while shear fractures form in even thicker layers. The absence of shear fractures in thinner layers may be explained by a geometry-related compressive stress decrease in the pinches. This results in localised shear evolving only in thicker layers. Complete layer thickness: fracture distribution, failure mode, and macroscopically resolved extension.


Introduction
Boudinage is the periodic failure -brittle or ductile -of competent layers in a mechanically layered material under coupled layer-parallel extension (e.g. Ramberg, 1955, Goscombe et al., 2004. For a discussion of definition of nomenclature see Urai et al. (2001) and Kenis et al. (2004Kenis et al. ( , 2005. Boudins provide a rich archive to understand the deformation history in the mechanically layered crust. First, boudins occur on all scales, from crustal-scale extension forming island chains (Jolivet et al., 2004), to km-scale boudins in thick layers (Strozyk et al., 2014), down to boudinage of brittle grains in ductile rocks. Second, their geometry depends on paleo-rheology, strain history, and initial geometry (Ramberg, 1955, Ramsay, 1968, Ramsay & Huber, 1983, Kenis et al., 2005. The importance of boudins for understanding deformation histories, as well as their geometric variability has been recognised since more than a hundred years (Ramsay, 1881, Lohest, 1909, Ramberg, 1955, Ramsay, 1968, Ramsay & Huber, 1983. A unifying taxonomy has been proposed by Goscombe et al. (2004), who subdivide boudins into five end members, in particular: ductile drawn boudins (also known as pinch-and-swell), brittle torn boudins separated by extensional (mode 1) fractures, bookshelf-rotated domino boudins separated by straight antithetic shear fractures, or gash boudins separated by sigmoidal or forked shear veins and fractures, and finally shearband boudins resulting from synthetic shearing. Our understanding of the mechanisms controlling boudinage and the resulting geometries has greatly improved in the past two decades, especially through analogue and numerical modelling. Analogue models (Zulauf & Zulauf, 2005, Zulauf et al., 2020a and numerical models , Virgo et al., 2013, Virgo et al., 2014, Virgo et al., 2016 reveal how complex boudin interference structures can evolve from multiphase and/or high-strain deformation. Numerical models have also shown to be a valuable tool in exploring the rheological control on boudin geometry and internal stress distribution (Bai et al., 2000, Schmalholz et al., 2008, Schöpfer et al., 2011, Abe & Urai, 2012, Moulas et al., 2014, Schöpfer et al., 2017. Nevertheless, we still do not fully understand the mechanics of boudin formation, especially the factors controlling boudin localisation. In this study we build on a high-resolution 3D dataset of an amphibolite layer that has experienced multiple phases of brittle and ductile boudinage . In this dataset, a sequence of boudinage phases has been identified in combination with extensive field studies . Furthermore, we showed that boudin geometry in 3D is much more complex than anticipated from 2D outcrop. Notably, the orientation and mode of inter-boudin fractures can change along strike, multiphase boudinage can lead to significant reactivation, strain recorded by boudinage is distributed heterogeneously, and the 3D dataset revealed an additional boudin set that was not obvious in the field . These studies show that, in line with previous studies (Schöpfer et al., 2011, Marques et al., 2012, Virgo et al., 2014, Zulauf et al., 2014, the role of structural inheritance and the effects of layer thickness on brittle boudin geometry and failure mode are still poorly understood. In this study we further investigate the effects of layer thickness on brittle boudinage. Specifically, we investigate (1) the influence of layer thickness on boudin width during polyphase deformation, (2) how layer thickness effects fracture mode during brittle boudinage, and (3) potential explanations for the correlation between layer thickness and resolved strain identified previously . Finally, we integrate the additional boudin generations discovered in the high-resolution dataset into the exhumation and cooling history of Naxos.

Boudinage Formation and Geometries
Brittle boudins evolve through sequential fracture infilling, i.e. fractures forming in between existing fractures, until the layer is saturated with fractures and a steady-state boudin width is reached, which is generally thought to depend on layer thickness, coupling between matrix and boudinaged layer, effective stress, as well as material properties of both layers, such as Young's modulus of layer and matrix, Poisson's ratio of layer and matrix, matrix viscosity, and layer tensile strength (Ramberg, 1955, Price, 1966, Mandal et al., 1994, Bai et al., 2000, Schöpfer et al., 2011, Marques et al., 2012. Geometrical complexity such as splay and fractures clustering in arrays (Olson, 2004, Abe & Urai, 2012 develops during the early stage of infilling, while maximum tensile stress is still heterogeneously distributed along the layer (Schöpfer et al., 2011). Once saturation is reached, fracture density should correlate directly with rheology and stress state in the competent layer (Bai & Pollard, 2000b, Schöpfer et al., 2011, Zulauf et al., 2020a. When coupling between matrix and competent layer is strong, interfacial slip is suppressed, and stress is redistributed inside the boudins with a critical aspect ratio of about 1. At this stage, further boudin segmentation is only possible through additional weakening factors such as material defects or elevated fluid pressure (Bai et al., 2000, Bai & Pollard, 2000b, Schöpfer et al., 2011.
In order to reconstruct the full deformation history recorded in boudin structures it is crucial to investigate their geometry in three dimensions, as only sectional and map views combined can reveal the full complexity of boudin structures (Zulauf & Zulauf, 2005, Zulauf et al., 2011b, Marques et al., 2012, Zulauf et al., 2020a. 3D boudin analysis has shown that identifying multistage deformation with the same failure mode is in general not possible from a 2D section alone . Similarly, tablet boudins (of which chocolate tablet boudins are a subtype) formed either in one episode of pure flattening (Ghosh, 1988, Zulauf et al., 2011b, in nonplane strain (Ghosh, 1988, or through polyphase deformation (i.e. the classical chocolate tablet boudins, Ghosh, 1988, Zulauf et al., 2011a, Zulauf et al., 2014)-these look identical in 2D sections but only reveal their deformation history in map view. Finite strain can also be underestimated when secondary reworking of boudins that have been parallelised to the instantaneous extension axis by in-plane folding goes unnoticed (Zulauf & Zulauf, 2005, Zulauf et al., 2020a.

Geological Background
The island of Naxos, in the centre of the Aegean Sea, offers a unique opportunity to study multiphase boudinage structures in three dimensions. Marbles, mined for their clear white colour and large crystal size, are intercalated with centimetre to decimetre thick amphibolite layers and pegmatite dykes, and have undergone multiphase deformation (Schenk et al., 2007. Naxos is a metamorphic core complex exhumed by a crustal-scale, N-dipping detachment with top-to-the-N sense of shear (Lister et al., 1984, Urai et al., 1990. Extensional tectonics in the back-arc basin are a manifestation of slab rollback following southward migration of early Alpine subduction in the Hellenic arc (Fytikas et al., 1984, Jolivet & Brun, 2010, Jolivet et al., 2013. The metamorphic core complex hosts a central migmatite dome subdivided in three sub-domes, with several marble enclaves and rafts in and around the high-strain zone separating the southern from the central sub-dome (Vanderhaeghe, 2004, Kruckenberg et al., 2011. Deformation structures in these marble rafts indicate multiphase E-W shortening evidenced by two generations of ductile drawn boudins isoclinal folding, and five generations of brittle boudins in the amphibolite layers . The pegmatites, which intruded the marble after ductile deformation in the amphibolite, only display two generations of brittle boudins (Schenk et al., 2007). The structures within the marble rafts are shown to post-date development of similarly coupled folding and boudinage in the metasediments surrounding the migmatite dome (Urai et al., 1990, Buick, 1991.

Sample & Methods
We present data from serial sectioning of a mined, wedge-shaped marble block with a height of 37 cm and face lengths of 250, 224, and 140 cm ( Fig. 1 in Von . It contains a single boudinaged amphibolite layer with variable thickness resulting from earlier drawn boudinage, cut by brittle boudins with axes subnormal to the drawn boudins. It was recovered from the marble mine south of Kinidaros on the island of Naxos, Greece (37°05'12.6" N, 25°28'20.0" E), where a large marble raft from the high-strain zone subdividing the migmatite dome is mined. The exact in-situ orientation of the block is unknown, but can be constrained by comparing the boudin structures to the field data from Virgo et al. (2018): the amphibolite layer had to be vertical and striking N-S. Amphibolite composition is dominated by hornblende and plagioclase, with biotite, muscovite, and minor titanite (Schenk et al., 2007).
We employed 'outcrop-scale tomography  to generate a 3D model of the boudinaged layer. The internal structure of the metre-sized sample is reconstructed from slices cut with an industrial diamond saw. The sample was cut into slabs perpendicular to the axes of brittle boudins, based on observations in the field. The slabs are about 2 cm thick, with a kerf loss of about 1 cm. Both faces of each slab were digitised using a Nikon D800 camera with a Nikkor Micro 60 mm f/2.8D lens, resulting in an image resolution of 0.2 mm/pixel and negligible distortion of 1 mm over a frame width of 1450 mm. The boudin structure was then reconstructed from the spatially referenced images in Midland Valley's (now Petroleum Experts) Move 2016 suite. In each slice (i.e. image), we mapped fault traces (classified according to failure mode), layer top and bottom, inter-boudin width (measuring dilatancy), and then correlated these across slices ( Figure 1). In total, 5087 fracture traces created a 3D network of 841 fractures planes. Details of the interpretation were tested by dissolving the marble in selected slabs using hydrochloric acid. Photogrammetry models of thus exposed amphibolite were integrated into the structural 3D model of the entire sample and confirmed our interpretation. The effects of thickness on fracture spacing, failure mode and resolved strain are quantitatively analysed in synthetic sections calculated from the 3D model and striking perpendicular to the thickness gradient of the amphibolite layer (4° clockwise misorientation to the slices, Figure 1). Microstructural analysis was performed using the Petroscan Virtual microscope (Virgo et al., 2016).

Results
In the following paragraphs, we show our observation in the reference frame of the sample, with the amphibolite layer being horizontal. Orientations in this reference frame are indicated by an asterisk (*). We note that the amphibolite layer was oriented sub-vertically with the E-W* axis trending N-S, as explained above.

Observed Structures
Amphibolite layer thickness varies across the sample due to early drawn boudinage. Drawn boudins are present in two generations , with the older, long wavelength (lλ) boudins responsible for the metre-scale thickness gradient from the long end to the short end of the sample. Their wavelength exceeds the sample width and is probably around 2-3 m, as observed in the field . Due to their long wavelength the shape of lλ boudins is not well constrained, but can be assumed to be cylindrical at the sample scale. They are overprinted by a second generation (sλ) of drawn boudins, which have a wavelength of ca. 10-20 cm ( Figure 1). The sλ boudins are not fully separated but form an anastomosing network of amphibolite interrupted by lenticular pinched-out areas of up to 1 m length and a few decimetre width (see Figures 9, 10, 14 in Virgo et al., 2018 for additional field images). Torn boudins with apertures up to 3 mm striking 178 ± 15°*, 417 inter-boudin fracture planes. They are the most abundant fractures in the block.
Veins: Chlorite filled veins striking 18 ± 14°*, frequency not recorded as they often do not form continuous planes in layer-normal direction. They are consistently protruding into the marble (i.e. thicker mechanical unit), have an aperture of 1-2 mm, a layer-normal length of 2-5 cm and a consistent spacing of 5-10 cm. In thick amphibolite they tend to be pinned to older, well-developed fractures.
Chlorite slickenfibres developed in many sheared inter-boudin zones (or fractures), recording sense of motion, dilatancy, and periodicity of deformation. Shear motion is pure dip-slip, with commonly a small extensional component (i.e. hybrid fractures). Reactivation as extensional fractures is frequent and evidenced by ruptured or new slickenfibres ( Figure 3). Both shear and extensional fractures are reactivated by chlorite veins in inter-boudin zones and layer-normal chlorite veins above and below shear fractures in the marble. The marble also has a lighter colour to either side of the amphibolite layer for as far as the chlorite veins protrude into it. The different boudin sets are strongly interconnected, with frequent bridging or reactivation and extension of older fractures by younger ones. Interaction between the NNW*, NNE* and N* sets is almost exclusively characterised by Y-nodes (i.e., younger boudins terminating at older boudins).

Effects of Thickness
The large-scale thickness variation of the amphibolite layer originating from early drawn boudinage provides an opportunity to investigate the effect of layer thickness on brittle boudinage. Other potential controls on boudin geometry like the amphibolite's mechanical properties, grain size, deformation history as well as external factors like far-field stress (and hence marble flow) can be assumed to remain constant within the sample for any given boudin generation, leaving amphibolite thickness as the only major variable. Fracture geometry and distribution reveal three parameters correlating with layer thickness: fracture distribution, failure mode, and macroscopically resolved extension.

Fracture Spacing
Boudin widths in the model (i.e. fracture spacing) can be analysed using two methods: A more intuitive approach is to calculate spacing from the number of fractures of each generation counted in sections that are perpendicular to the thickness gradient of the amphibolite layer (referred to as top-down approach). This has the advantage that we can calculate spacing for each boudin generation and thus observe the influence of inheritance. However, thickness-normal sections are synthetic and only contain intersections with the interpolated and classified fracture planes (grouped into generations), potentially introducing bias. These drawbacks can be circumvented by the second approach, where we directly count the mapped 2D fracture traces in the original slabs, disregarding interpolation across slabs and classification (referred to as bottomup approach). This, however, oversimplifies the system as it does not account for finite strain rotation between different boudin generations, forcing new localisation instead of reactivation of older fractures. Also, the approach might overestimate fracture density where fractures interact or splay. Finally, since 3D fracture orientation is unknown and changes from fracture to fracture, we cannot correct the spacing calculated in the bottom-up approach for the length-overestimation resulting from measuring on non-perpendicular sections. The resulting overestimation of spacing is significant only for the low-density NW* set (88.7 %) but acceptable for all other boudin sets (NNW*: 8.6 %, NNE*: 0.4 %, N*: 0.1 %, Veins: 5.1 %).
Results from both approaches are very similar for extension fractures but less so for shear fractures (Figure 4). NNE* shear fracture spacing (top-down) follows a negative power-law correlation with thickness, with an exponent of -0.8 ± 0.11 and an 2 of 0.6. Conversely, shear fracture spacing in the bottom-up approach shows a positive correlation with thickness, like N* fractures in the top-down approach, as well as extensional and all mapped fractures in the bottom-up approach. Fitting parameters are summarised in Table  1. Extensional fracture spacing (N* in top-down and extensional in bottom-up) are potentially uncorrelated with layer thickness as 2 for both fits are < 0.1. The average aspect ratio of boudins (fracture spacing over layer thickness) is quite low, at 0.3 ± 0.11. This matches the qualitative observation that, at least where the amphibolite is thick, wider boudins with an aspect ratio of about 1 are separated by clusters of several fractures that are sometimes interconnected (Figure 1). The high standard deviation of fracture spacing measured in subsets of fractures (i.e. NNE*, N*, shear, or extensional fractures) is reduced if data for shear and extensional fractures are summed up to calculate overall fracture spacing (9 mm & 16 mm vs 4 mm, respectively). This indicates that there has been reactivation. To test this, we plotted fracture spacing and amphibolite thickness along the thickness gradient ( Figure 4C). The plot reveals a striking negative spatial correlation between the two dominant boudin generations in the sample, N* (extensional) and NNE* (shear) boudins. N* fractures preferentially develop in areas with less pre-existing (NNE*) fractures. The trends of both boudin generations also shows a distinct sinuosity, not related to any measurable thickness variation, which is strongly reduced if their spacing is averaged.

Failure Mode
Extension in the amphibolite layer is dominantly accommodated by slip on shear fractures where the layer is thick, as opposed to dilatancy of extensional fractures where it is thin ( Figure 5A). This observation is corroborated by the dip vs thickness ( Figure 5B) and extension vs thickness data ( Figure 5C), both showing increasing shear dominance with increasing layer thickness. Dip angles of NNE* and NNW* inter-boudin fractures, which both have a distinct shear component, display a very similar negative correlation with thickness with an average factor of about -0.35°/mm. Dip angles of the purely extensional N* and NW* fractures remain largely constant over the entire thickness range. Dip angles are used as proxies for the failure mode, with the underlying assumption that extensional fractures should develop at ca. 90° to the layer interface and this angle decreases as the shear component of deformation increases.

Brittle Extension
Extension localised in brittle boudinage is either accommodated by purely layer parallel dilatancy or by heave, the layer parallel component of shear slip. Slickenfibres in the etched blocks confirm the movement between boudins to be purely dip-slip. Dilatancy is directly mapped during section interpretation and heave is the difference between finite boudin train length, boudin exterior, and dilatancy (see Von . Both are normalised to finite boudin train length and provide a direct approximation of finite strain accommodated by brittle boudinage ( Figure 5C). While dilatancy is essentially thickness invariant, heave correlates with thickness. As a result, their sum, total brittle extension also correlates with layer thickness, ranging from ca. 5 % extension at 5 mm thickness to 16 % at 32 mm. The power-law fit for heave suggests a lower boundary for the development of shear fractures around 10 mm thickness. Further, dilatancy below this thickness is exclusively accommodated by chlorite veins, which are controlled by a different mechanical unit (they extend into the marble) than the previous boudin generations. Hence, we can consider 10 mm to be the threshold below which no layer-bound fractures localise in our sample, and indeed, unless linked to a chlorite vein, we can observe deformation progressively delocalising below 10 mm layer thickness. Just below this threshold we can still observe some necks and shear zones at angles of ca. 35-45° to the layer, which disappear as layer thickness approaches the grain scale ( Figure 6A). The correlation between fracturing and thickness is also reflected in failure mode and fracture spacing, notably the negative power-law trend of the NNE* set ( Figure 4A).

Microstructure
Thin sections from thick and thin amphibolite show similar grain size distribution (plagioclase: 0.2-0.5 mm, amphibole: 0.3-1 mm), degree of crystal plastic deformation, frequency of deformation twins in plagioclase and calcite, and frequency of grain cracking. We did not measure grain size distribution (e.g., by EBSD) nor study crystal plastic deformation in the surrounding marble (this could be done using the mechanical twins which are ubiquitous): the samples are available for further study. Only two differences are apparent: The concentration of calcite and microscopic shear planes. While calcite is absent in thick parts of the amphibolite, small, disseminated pockets (which could well be calcite-filled dilatant microfractures) of up to few hundred micrometres across become more common where the layer thins. They have highly angular outlines with wall rock minerals commonly protruding into them, hinting at their secondary origin.
Distributed microscopic strain is accommodated by shear planes across the entire amphibolite layer ( Figure  6B-E). They are characterised by consistent alignment of grain boundaries along a single surface, concentration of small, rounded grains, biotite, and chlorite along the aligned boundaries, as well as occasional intragranular deformation, either as slip along cleavage planes, asperity removal, or crystal plastic deformation. The latter is most commonly observed as undulose extinction and subgrain boundaries in calcite grains of the marble, but can sometimes also afflict plagioclase grains. Our observations do not reveal if the grain boundary alignment is an effect of asperity abrasion on the micrometre scale or crystal-plastic processes. When abutting against larger grains, shear planes can bifurcate around the grain and reunite, or splay and continue as multiple planes. Less frequently in the amphibolite, but commonly in the coarser marble, they penetrate and terminate within a grain. Where identifiable, displacement on these planes does not exceed a few micrometres. These shear planes are abundant throughout the whole layer in thick amphibolite, irrespective of proximity to fractures. In thinner amphibolite they are only abundant in and around macroscopically identifiable sites of deformation-either more localised shear bands or less localised necks. As such, they are another manifestation of brittle strain decreasing with layer thickness.

Spacing
Due to the thickness variation from old drawn boudins, our sample provides a unique opportunity to evaluate the effects of layer thickness on boudinage under otherwise constant conditions of temperature, pressure and far field strain rate. Results show a clear relationship between thickness and failure mode ( Figure  5) as well as total brittle extension as measured by fracture displacements ( Figure 5C). However, the relationship between layer thickness and fracture spacing is far more ambiguous as would be expected from previous studies (e.g. Mandal et al., 1994, Bai et al., 2000, Schöpfer et al., 2011, Marques et al., 2012. Fracture spacing in the individual boudin generations shows only weak correlation with thickness in the case of the NNE* set, but not the N* set ( Figure 4A). Moreover, spacing of NNE* fractures correlates negatively with thickness, which implies that strain is not distributed homogeneously through the boudinaged layer. The negative trend is, however, consistent with resolved strain increasing with layer thickness. The high variance in spacing of both NNE* and N* sets, and thus the weak correlation, can partly be explained by complex fracture geometries with abundant splaying and branching, as well as reactivation of pre-existing fractures being favoured in the already saturated parts of the layer. Reactivation can be expected for the low misorientation between the NNE* and the N* set, especially given the potentially high rheological contrast between incompetent marble and competent amphibolite (Virgo et al., 2014). Indeed, interaction between the NNW*, NNE* and N* sets is almost exclusively characterised by Y-nodes, suggesting that newly developing fractures readily propagated towards and merged with existing fractures. The significance of structural inheritance becomes evident when plotting spacing along a thickness section through the amphibolite ( Figure  4C), showing the alternating dominance of either the older NNE* or the younger N* set, while the average combined spacing stays relatively constant. This highlights how reactivation is an effective way to suppress further fracture nucleation where fracture density is already high (Passchier et al., 2021, Prabhakaran et al., 2021a, Prabhakaran et al., 2021b. The NNE* fractures were reactivated during the subsequent N* deformation, and new fractures only developed where the former were scarce. The considerably lower standard deviation of fracture spacing across all boudin generations ( Figure 4B) compared to individual boudin generations ( Figure 4A) further supports this hypothesis. Extensional fracture spacing shows no significant correlation with thickness, and spacing between shear fractures as well as between all fractures, not distinguished between fracture modes or boudin generation, follows a low-exponent power law or linear trend. Assuming a fracture infill model (Bai et al., 2000, Bai & Pollard, 2000a, Bai & Pollard, 2000b, new fracture nucleation should stop as boudin aspect ratio approaches a critical value (usually around 1), and additional strain should be accommodated in the inter-boudin zones (i.e., by the fractures in our case). However, the measured boudin aspect ratio of only about 0.3 is significantly lower than the expected critical aspect ratio of ca. 1 (Bai et al., 2000), mainly due to the clustering of fractures and splays. Our boudin sections resemble the model results with high strain and interfacial friction coefficient in Fig. 4 of Schöpfer et al. (2011). Clusters of fractures can form during the infill phase of fracturing when tensile stress is not yet bound to the boudin centre, and flaws in the direct proximity of an existing fracture can propagate through the layer or link to the existing fracture (Bai & Pollard, 2000a, Schöpfer et al., 2011. Our data support an overall weak correlation between fracture spacing and layer thickness, but other factors such as structural inheritance and geometrical complexity arising from fracture nucleation and propagation are more important for controlling fracture distribution during individual generations of deformation.

Failure Mode
The variation in strain accommodated by ductile and brittle boudinage can be explained by differential embrittlement and viscosity increase during exhumation and cooling, as well as, in the case of brittle boudins, the geometrical constraints on failure mode arising from earlier drawn boudinage. Boudin size decreases in each consecutive boudin generation (i.e., from 24 inter-boudin fractures with negligible aperture in the early NW* generation to 178 inter-boudin fractures with ≤3 mm aperture in the late N* generation, Figure 2). Progressive early embrittlement of the amphibolite is recorded by the decreasing boudin size from the lλ (2-3 m wavelength) to sλ (10-20 cm wavelength) generation and ultimately the transition from ductile to brittle boudinage. Continued viscosity increase of the marble then lead to increasing fracture density in each subsequent brittle boudin generation as the viscosity contrast between matrix and layer decreased and interfacial coupling increased (Schöpfer et al., 2011, Zulauf et al., 2020a. Finally, the marble was embrittled at least locally when it was fractured by chlorite veins. This confirms the continuous trend of cooling and exhumation proposed in Virgo et al. (2018). As a consequence, earlier localisation was more likely to develop as shear fractures (NNW* and NNE* sets) under the higher confining pressure, while the later N* and vein sets developed as extensional opening mode fractures (e.g. Schöpfer et al., 2007, Abe & Urai, 2012, Kettermann & Urai, 2015. This could be explained by the increased viscosity of the marble limiting material flow around fractures (Schueller et al., 2005), preventing significant block rotation from shear fracturing. A reduced confining stress would also explain why the chlorite veins protrude into the marble, as they could propagate into the then relatively competent marble after nucleating in the brittle amphibolite. The precipitation of chlorite veins also evidences the presence of fluids, although there is not much evidence for retrograde reactions in the amphibolite itself. However, a detailed petrographic study was not attempted in this project. A late-stage fluid pressure pulse reducing effective stress in the amphibolite and surrounding marble could have promoted fracture propagation into the marble (see also Schenk et al., 2007). Further, shear fractures in our sample preferentially localised in drawn boudins (or swells) as can be seen in Figure 5. During the earlier NNW* boudin generation, shear fractures only developed in the northern* swell with a fairly circular crosssection (low aspect ratio), while during the later NNE* boudin generation shear fractures also developed in the more drawn out southern* swell (high aspect ratio) (Figure 2 & Figure 5). Mandal et al. (2001) show how tensile stress in brittle inclusions increases linearly with their aspect ratio, while compressive stress decreases, making low aspect ratio swells more likely to fail in shear while the high aspect ratio pinches fail in tension. Hence, the geometric constraints from early drawn boudinage dictated stress partitioning in later phases of brittle deformation. Shear fractures first nucleated in the swells with the lowest aspect ratio and then progressively propagated into the thinner and more drawn-out parts of the amphibolite layer.

Strain
The asymptotic increase of fracture spacing for NNE* shear fractures in thinner amphibolite ( Figure 4A) is emphasising how thinner layers could not accommodate shear fractures. However, the 'missing' strain was seemingly not transferred to extensional fractures either. The same trend is apparent in the rotation of fracture dip in the shear dominated NNW* and NNE* sets towards layer normal extensional fractures in thinner amphibolite, while the tensile-dominated N* and NW* sets are thickness invariant ( Figure 5B). The gradual change in orientation is distinct from the sudden changes in orientation and failure mode resulting from lateral fracture propagation by reactivation. The correlation between brittly resolved strain and thickness is further corroborated by measured macroscopic extension ( Figure 5C). Correlation of brittle extension with thickness can mainly be ascribed to the behaviour of shear fractures ( Figure 5C). However, we note that the marble has not been stretched homogenously, and we can determine local strain variations on the order of 10 % at the sample scale (ca. 4.5 m 2 ). Since below 10 mm layer thickness only chlorite veins disrupt the layer, we infer that layer-bound brittle deformation was replaced by macroscopically distributed deformation (Figure 7). This transition to macroscopically distributed deformation in thin layers can also be seen in field data from Naxos, where adjacent layers of different thickness have very different fracture densities and associated macroscopic brittle extension. Fig. 7 in Virgo et al. (2018) and Fig. 4 in Schenk et al. (2007) show examples of thin amphibolite layers following the stepped boudin exterior of thicker neighbouring amphibolite layers or pegmatite dykes. These thin 'sidekick' layers must accommodate even more strain than their thicker neighbours but show no macroscopically localised deformation. Consequently, these layers must have accommodated extension by distributed processes as overall fracture density decreases (i.e. spacing increases). Below a thickness of ca. 10 mm, brittle deformation is limited to the chlorite veins, which are however not constrained to the amphibolite layer and were most likely controlled by a thicker mechanical unit including the lighter coloured marble directly adjacent to the amphibolite (Figure 1). Just below the thickness threshold for brittle localisation, deformation could still 'localise' in cm-wide zones that manifest macroscopically by necking, but without separation. Microscopically, these zones contain abundant shear planes (Figure 6), mostly following and aligning grain boundaries. However, these surfaces are much more abundant in thicker amphibolite, where they are not limited to the vicinity of fractures but are also common in boudin centres ( Figure 6E). Also, the deformation mechanism governing these shear planes remains elusive. They show no evidence for cataclasis, compaction, or mineral alteration. The distribution of the shear planes (clustering in necks from brittle boudinage) suggests that distributed brittle strain is also proportional to layer thickness, which means that these shear planes are not accommodating 'missing' extension in the pinches. Increased concentration of calcite in thin amphibolite could indicate intergranular fracturing and subsequent granular flow as a possible mechanism for distributed deformation, as described for the drawn boudins by Virgo et al. (2018). However, to preserve the high grain angularity seen in the amphibolite, considerable dilatation would have been necessary, which would require elevated fluid pressure that should have manifested in much more significant mineralisation (Stel, 1981). The absence of layer-bound fractures in very thin amphibolite is possibly an effect of intergranular fracturing (Abe & Urai, 2012, Spruženiece et al., 2021 saturating the layer before cracks could coalesce into localised fractures. This is not surprising considering that the layer is thinner than the expected shear band thickness of 10 -20 grains across (e.g. Mühlhaus & Vardoulakis, 1987).

Conclusions
The Naxos amphibolite boudins offer exciting insights into the evolution of a polyphase boudin system under progressively decreasing temperature and pressure. Our data reveal three unexpected effects: (1) locally, boudin width is controlled by inheritance rather than layer thickness, and (2) fracture mode as well as (3) brittly resolved extension are proportional to layer thickness. These observations are obvious when analysing boudins in 3D, but would probably remain unnoticed in a description of a typical 2D outcrop in the field. Thus, our findings also highlight the necessity for characterising such complex systems in 3D, leading to the following geometric and mechanical findings, specifically: 1. Fracture density and failure mode record the progressive cooling and embrittlement of marble and amphibolite during exhumation. 2. Failure mode of brittle boudins varies with layer thickness. Shear fractures were only accommodated in swells with a low aspect ratio cross-section, while extensional fractures could also localise in pinches.
3. Brittle strain delocalised with decreasing layer thickness. Below the threshold of 10-20 grains thickness, localised fracturing is absent and microscopic shear planes become constrained to shear zones and necks. 4. As a consequence of delocalisation, resolved strain in the boudinaged layer correlates with thickness and thin domains record less strain than thick domains. No process accommodating this 'missing' brittle strain could be identified. Presuming the layer stretched homogeneously, the accommodating process remains elusive. 5. Final total fracture spacing (or boudin width) is controlled by layer thickness. However, within individual boudin generations, influences such as structural inheritance, geometric complexity from fracture nucleation, and the effects of failure mode variation are dominant.