Mixed Brittle and Viscous Strain Localization in Pelagic Sediments Seaward of the Hikurangi Margin, New Zealand

Calcareous‐pelagic input sediments are present at several subduction zones and deform differently to their siliciclastic counterparts. We investigate deformation in calcareous‐pelagic sediments drilled ∼20 km seaward of the Hikurangi megathrust toe at Site U1520 during International Ocean Discovery Program (IODP) Expeditions 372 and 375. Clusters of normal faults and subhorizontal stylolites in the sediments indicate both brittle faulting and viscous pressure solution operated at <850 m below sea floor. Stylolite frequency and vertical shortening estimated using stylolite mass loss, porosity change, and distribution increase with carbonate content. We then use U1520 borehole data to constrain a P‐T‐t history for the sediments and apply an experimentally derived pressure solution model to compare with strains calculated from stylolites. Modeled strains fail to replicate stylolite‐hosted strain distribution or magnitude, but comparison shows porosity, composition, and grain‐scale effects in diffusivity and mass transfer pathway width likely exert a strong influence on pressure solution localization and strain rate. Stylolite and fault clusters concentrate clay in these sediments, creating weak volumes of clay within carbonates, that may localize slip where the plate interface intersects the carbonates at <5‐km depth. Plate interface slip character and rheology will be influenced by the deformation of intermixed phyllosilicates and calcite, occurring by variably stable frictional slip and pressure solution of calcite. Pressure solution of calcite is therefore important at the shallow plate interface, waning at the base of the slow‐slipping zone because calcite solubility is low at temperatures >150°C where frictional (possibly seismic) slip likely predominates.


Introduction
The mechanical behavior of carbonate-rich sediments during subduction is not well understood, despite their recognition at several margins around the world (Moore & Mascle, 1990;Morris et al., 2006;Wallace et al., 2019). Rock deformation experiments have, however, shown that carbonates deform very differently to siliciclastic sand and mudstones typically studied on subduction margins (Boulton et al., 2019;Ikari et al., 2013;Kurzawski et al., 2018;Rabinowitz et al., 2018). For pressure (P) and temperature (T) conditions at <15-km depth (T < 200°C) in subduction margins, siliciclastic lithologies deform dominantly by cataclasis (combinations of fracture and frictional sliding) with slower pressure solution, whereas carbonates undergo more appreciable deformation by pressure solution and crystal-plastic deformation, possibly associated with seismic slip (Gratier et al., 2013;Kennedy & White, 2001;Verberne et al., 2013). The activity of pressure solution at low temperature in carbonates could have significant implications for the rheological behavior of subduction thrust interfaces that contain calcareous sediments.
Several factors influence the rate of pressure solution, including clay content (Aharonov & Katsman, 2009;Hickman & Evans, 1995), grain surface area (Gratier et al., 1999), mineral solubility (Rutter, 1976), grain contact area and microstructure (including variations due to porosity and grain boundary healing) (Croizé et al., 2010;Ebner et al., 2010;van den Ende et al., 2019), fluid flow and composition (Lehner, 1995;Zhang & Spiers, 2005), and gradients in normal stress (Rutter, 1976). The constitutive equations for creep rate by pressure solution are distinct when rate limited by either dissolution, diffusion, or precipitation, and each is consequently described by a different model (Gratier et al., 2013;Raj, 1982;Rutter, 1976). Pressure solution on stylolites at shallow depth represents viscous deformation in a depth range typically associated with brittle fracture and frictional sliding (<15 km Paterson & Wong, 2005;Tada & Siever, 1989). At constant environmental conditions, the rate of pressure solution is dominantly controlled by mineral solubility and grain size, facilitating greater ductility of the upper crust in lithologies containing more soluble minerals (Gratier et al., 1999;Renard et al., 2000). The brittle and viscous components of deformation in the upper crust therefore vary in relative intensity with lithology as well as depth and temperature. Where elevated rates of pressure solution form stylolites, changes in properties such as clay content and porosity affect localization and frictional slip properties of faults (Baud et al., 2016;Tondi et al., 2006;Viti et al., 2014;Watkinson & Ward, 2006).
International Ocean Discovery Program (IODP) Expeditions 372 and 375 drilled and sampled five sites on a transect across the Hikurangi Margin, New Zealand, between November 2017 and May 2018 (Wallace et al., 2019). We study calcareous-pelagic chalks and marls in the incoming sedimentary sequence sampled at Site U5120. We measure fault and stylolite distribution and make microstructural observations to determine the extent and character of pre-subduction deformation before calculating uniaxial strain from pressure solution using chemical mass loss estimates and stylolite frequency. A possible relationship between CaCO 3 content and stylolite frequency could provide an insight into a key natural control on stylolite development.
The abundant data gathered during IODP Expedition 375 provide a unique opportunity to constrain the P-Tt history for the pelagic sequence at Site U1520, which we use to apply published models of intergranular pressure solution (Rutter, 1976(Rutter, , 1983. We then compare these results to those calculated from stylolite frequency earlier in this study and find agreement requires sustained high fluid pressures or larger length scales of material transport than the average grain size of lithologies hosting stylolites. Finally, we extend this model to the future subduction of the sediments and discuss the rheological impact of ongoing pressure solution and pre-subduction strain weakening on current subduction. on the northern Hikurangi Margin accommodate up to 250 mm of slip over several days or weeks and recur every 2-5 years Wallace et al., 2016). As well as slow slip, the northern part of the Hikurangi Margin is characterized by episodic seismicity of variable character, including large tsunami earthquakes (≥M w 7.2, Clark et al., 2019;Doser & Webb, 2003;Wallace et al., 2014), microseismicity (Delahaye et al., 2009), and tectonic tremor (Todd & Schwartz, 2016;Yabe et al., 2014).
Seismicity at the Hikurangi Margin is distributed throughout the sedimentary upper plate and oceanic lower plate, with tsunamigenic earthquakes occurring on the plate interface (Clark et al., 2019;Shaddox & Schwartz, 2019;Todd et al., 2018). The volume hosting slow slip likely encompasses the plate interface and shallow earthquake hypocenters, extending from near the trench to 15-to 20-km depth (Wallace et al., 2016), possibly hosted in heterogeneous lithologies characterized on the incoming plate (Barnes et al., 2020;Wallace et al., 2019). Seismic transects across the margin have shown relatively thin sediment cover on the incoming plate (≤1 km) between large seamounts (Barker et al., 2009;Davy et al., 2008). High-amplitude reflections are present in both the incoming plate and downdip of the slow slip-hosting volume,  Wallace et al. (2016) east of North Island, New Zealand (inset map). The red point is IODP Expedition 375 drill Site U1520, on the white line showing the 05CM-04 seismic section drawn below (c). The Hikurangi Trench is in red. (b) Lithostratigraphic log shows core from site U1520. Colors of units (yellow, green, brown, and blue), and contacts between units (orange, green, and brown) are the same on section (c) below. Schematic section (c) is redrawn from Barnes et al. (2020) and shows the relationship of the units cored at U1520 to the plate interface. Also shown are the décollement (thick black line), active faults (red lines), inactive faults (dark red lines), normal faults (thin black lines), base of volcanic basement (light blue line), approximate slow slipping zone, and the high-reflectivity zone (HRZ) as characterized by Bell et al. (2010).

10.1029/2019TC005965
Tectonics suggesting incoming material may be underthrust (Bell et al., 2010). The thin sediment cover, moderate convergence rate, rough incoming topography, and steep taper angle at this margin suggest it may be undergoing frontal erosion and is consequently likely to underthrust sediment rather than accrete it (Barker et al., 2009;Fagereng, 2011b).

IODP Drilling at Site U1520
Site U1520 is on the flank of the Tū ranganui Knoll seamount and therefore represents a slightly condensed incoming sedimentary sequence. Logging-while-drilling (LWD) data were collected at Site U1520 during Expedition 372 (Wallace et al., 2019). The same site was subsequently cored and wireline logged during Expedition 375, from 0 to 642.3 m below sea floor (mbsf) in hole U1520D, and from 646.0 to 1054.1 mbsf in hole U1520C (Wallace et al., 2019). Cores revealed a ∼1-km-thick sequence of heterogeneous sediments including turbiditic silts and muds (Units I-III), pelagic muds and carbonates (Unit IV), and volcaniclastic conglomerates with minor mudstones (Units V and VI; Figure 1b). In this work, we are concerned with the pelagic carbonates and carbonaceous muds which comprise Unit IV. Unit IV was found at depths of 509.82-848.45 mbsf (Wallace et al., 2019).
Site U1520 is approximately 20 km seaward of the toe of the Hikurangi megathrust ( Figure 1a). Assuming a convergence rate of ∼40 mm year −1 , the sediments intersected at Site U1520 will reach the toe of the thrust in approximately 0.5 Myr. Comparing the depth of Unit IV at Site U1520 with interpretations of active seismic data along line 05CM-04 , the shallower and deeper contacts likely correspond to the bases of seismic Units 5 and 8 (Barnes et al., 2020). These same horizons can be traced landward until truncated by the megathrust at depths of ∼4.5 and ∼5 km (Figure 1c), where the mud-and chalk-dominated rocks are inferred as a key protolith interval in which the plate interface will likely locate (Barnes et al., 2020).
Measurements of many physical and chemical properties were collected by LWD and shipboard core logging, details of which are described by Wallace et al. (2019). Here, we use data concerning the density, temperature, age, porosity, mineralogy, and lithology of sediments in Unit IV of Site U1520.
Relative mineral abundance from X-ray diffraction (XRD) is normalized using coulometric CaCO 3 content analysis (Wallace et al., 2019). Relative abundances of minerals determined by XRD were calculated for a simplified suite of four mineral groups: calcite, quartz, feldspars, and clays. Clays are measured as a bulk single mineral, though dominant species present in Unit IV are smectite and illite (see XRD methods in Wallace et al., 2019).

Conditions of Pre-Subduction Compaction on Stylolites
The physical properties of these sediments were logged by LWD during Expedition 372 and on extracted core during Expedition 375. From shipboard density measurements (i.e., moisture and density in Wallace et al., 2019), in situ effective vertical stress (σ ef f v ) was calculated as where there are n intervals of bulk density measurement spacing dz with measured bulk density ρ b . Gravitational acceleration (g) is equal to 9.81 ms −1 , and average seawater density (ρ w ) is assumed to be 1,024 kg m −3 (Nayar et al., 2016). Typical values for dz (i) are ∼1 m, though this varies with gaps in coring and core recovery. Hydrostatic (and, for comparison, near-lithostatic) pore fluid pressure is assumed throughout the studied depths. The results of the hydrostatic pore fluid pressure calculation (Figures 2a and 2b) show that effective stress in Unit IV increases approximately linearly from ∼4 MPa near the top (509.8 mbsf) to ∼8 MPa near the base of the unit (848.5 mbsf).
Five of the seven downhole temperature measurements conducted in the shallowest 250 m of hole U1520D form a linear trend corresponding to 38°C km −1 (Figure 2c). The other two values appear as outliers and are therefore ignored. Extrapolating a constant trend from these measurements to depths greater than the shallowest 250 mbsf, temperature increases from 19.4°C at 509.8 mbsf to 32.2°C at 848.5 mbsf (Figure 2c). A

10.1029/2019TC005965
Tectonics detailed discussion of geothermal gradient in this type of sediment is beyond the scope of this work, so we will focus on the conditions under which strain was likely accommodated in Unit IV. To constrain the depositional age of the sediments at Site U1520, we use an approximate age model based on nannofossil, planktonic foraminifera, and palaeomagnetic reversal ages (Figure 2d Wallace et al., 2019).

Characterising Stylolite and Fault Frequency and Texture
To determine pre-subduction mechanisms of sediment deformation, we counted faults and stylolites in core samples of lithostratigraphic Unit IV. Image scans of drill core (collected shipboard using the Section Half Imaging Logger, see core handling and analysis methods in Wallace et al., 2019) were used to visually identify stylolites and faults (e.g., Figure 3a). For the purposes of counting, stylolites were defined as typically dark, slightly wavy surfaces where material loss was apparent and concentration of insoluble materials was inferred. Faults were defined as closed or recrystallized fractures with visible apparent offset. The use of these definitions aided in discounting fine clay-rich sedimentary layers and drilling damage from stylolite and fault populations. Horizons of increased deformation measured here coincide with those noted by shipboard structural geologists based on direct observation of the core (Wallace et al., 2019). The depth where a fault or stylolite occurred within the core was recorded as the middle of an interval defined by the intersection between the planar structure and the two sides of the core liner. As noted by shipboard scientists measuring selected structures in 3D (Wallace et al., 2019), most stylolites are approximately horizontal.
Six samples were collected from throughout the pelagic carbonate sequence (see repository files for data collected from all samples) and set in epoxy resin to avoid breakage during thin section preparation. By careful cutting and grinding using silicon carbide grit slurry and a Logitech PM5 polisher with 0.3 μm aluminum oxide, polished thin sections were made from samples cut either perpendicular to bedding or parallel to lineation and perpendicular to fault surfaces. Polished thin sections were used to produce optical photomicrographs, backscatter electron (BSE) images and energy dispersive spectra (EDS). BSE and EDS data were captured using a Zeiss Sigma HD Field Emission Gun Analytical scanning electron microscope (SEM) at the School of Earth and Ocean Sciences in Cardiff University. BSE images were collected at 20 keV with an aperture of 120 μm and a working distance of 8.9 mm. We use the number of dark pixels per row of BSE image as a proxy for porosity as BSE returns no result when analyzing a void. Following Verberne and Spiers (2017) and Heilbronner and Barrett (2014), images were cropped, rotated, histogram matched using the "Match

Tectonics
Color" function in Adobe Photoshop, and Gaussian filtered with a Gaussian width of 4 before the algorithm of Otsu (1979) was applied to determine a darkness threshold between 0 and 255 corresponding to porosity. Pixels values below this value (52 for sample 20R4W and 50 for sample 19R1W) were counted, and their ratio to the total number of pixels per row was averaged over approximate areas to show a bulk porosity change.
Where cracks were present, we cropped the averaging area; both cropped and uncropped averaging areas within the stylolite are shown in Figure 4b. For EDS mapping of faults and stylolites, spot size was either 1.59, 2.49, or 2.66 μm, accelerating voltage was 20 keV, dwell time was 100 μs, and aperture size was 120 μm. For determination of element concentrations from SEM EDS images, relative weight percentages were calculated in Oxford Instruments AZTEC software. The element concentrations calculated are relative within a single mapped area.

Description of Stylolites and Faults
Macroscopically, more well-developed stylolites appear as composite structures comprising discrete dark seams, frequently linked with more subtle seams within the adjacent few centimeters of core. Some stylolitic surfaces within Unit IV have relatively steep apparent dips (>45°) and host minor shear offset on undulating dark surfaces. These possibly represent poorly developed early fault surfaces which have later undergone pressure solution or "healing" (Figure 3a). At the microscale, most seams are composed of several finer seams (Figures 3b and 3c), with their most intense shared anastomosing branch representing the macroscopically visible seam. Seams within these finer networks are relatively thin (<0.1 mm), undulating, and subparallel to bedding (Figure 3b). Calcareous microfossils of uniform size are common throughout Element maps show a Ca-rich protolith matrix with stylolite seams relatively enriched in magnesium, silicon, iron, and aluminum (Figures 3f and 3g). The bulk stylolite hosts increased silicon and magnesium while individual seams appear to dominantly concentrate iron (Figures 3f and 3g). Individual seams are of variable thickness and anastomose throughout the width of the stylolite to form a weakly connected network with scattered parallel fractures (Figures 3f and 3g). Fractures are open but host some aluminum, probably representing clays seen within seams in BSE images (Figures 3d and 3e). Open fractures could be a result of sample preparation, where clay was removed by polishing, or core extension, where extension parallel to the core axis occurs when retrieving core from depth. Because we interpret them as formed either by sample preparation or by drilling-induced damage, these open fractures are not included in any discussion or quantification of subseafloor deformation.
From BSE images and EDS maps, it is qualitatively clear that stylolites have lower porosity than their surroundings (Figures 3d-3g); mean measured porosities adjacent to two stylolites are between 0.41 and 0.77 ( Figure 4). Shipboard porosity measurements on bulk intact samples from within 1-m (19R1W) and 7-m (20R4W) depth of the samples shown in Figure 4 were between 0.36 and 0.42, below our estimates from image analysis ( Figure 4). With further dissolution to develop a stylolite further, porosity loss is likely to increase (Toussaint et al., 2018), suggesting that 0.4 is a more appropriate figure for mean porosity loss in moderately well-developed stylolites. This porosity loss is a mean change over anastomosing seams and interseam areas of the bulk stylolite, commonly hosting higher porosity (0.4-0.8) where undissolved shells are present and lower porosity (<0.2) where clays are concentrated ( Figure 4). Because pore-space within intact shells becomes connected to other porosity as soon as shells breakdown either mechanically or by dissolution, we include this fossil-hosted porosity in our measurements.  (Figures 5b-5d). Clay-rich material is present as relatively thick (≥3 mm) seams in the immediate vicinity of some slip surfaces (Figures 5a-5d). In SEM BSE images, it is clear that thicker clay-rich seams adjacent to faults are connected to similar, thinner, seams which run along fractures adjacent to the fault surface ( Figure 5e). EDS maps show enrichment of magnesium, iron, silicon, and aluminum throughout the clay-rich network, including the material comprising the fault surface ( Figure 5f).
Faults and stylolites rarely occur on the same surface within Unit IV, except where early faults have been exploited by pressure solution (Figure 3a). Where faults and stylolites occur within the same depth range, stylolites are typically intersected and offset by faults. Some faults host clay-rich material immediately

10.1029/2019TC005965
Tectonics adjacent to their slip surface, possibly indicating the existence of stylolites prior to faulting. The existence of stylolites adjacent to fault slip surfaces is uncertain as faults obscure the stylolite texture. These possible stylolites have therefore been excluded from stylolite frequency counts.

Downhole Frequency of Stylolites and Faults
Following the evidence for pre-subduction deformation on faults and stylolites, we now quantify their frequency and distribution within core samples. To analyze the distribution of structures as a function of lithology, we use a simplified lithostratigraphic sequence, ignoring minor lithologic intervals (<5 m thick) to get seven lithological units ( Figure 6).
Stylolites in Unit IV occur in distinct clusters within similar grain size sediments, increasing in intensity with depth. Clusters are more frequent and intense where CaCO 3 content is >50 wt%, most notably in chalk near the deeper contact with the volcaniclastic conglomerate ( Figure 6). The deepest chalk unit sits between 797.99 and 848.45 mbsf, contains >90-wt% CaCO 3 , and hosts ∼75% of stylolites in the section. Other significant peaks in stylolite frequency are in the two marl units (655.15-720.93 and 738.68-773.26 mbsf). Minor peaks in stylolite frequency occur throughout the deeper part of the sequence in the marl, chalk, and debris flow units ( Figure 6). Stylolites are very rare in the calcareous mudstone and conspicuously absent in the mudstone unit, which contains less CaCO 3 (<50 wt%). Similar to stylolites, faults are also clustered, but fault clusters are more distributed throughout the core than stylolite clusters.

Clustering Analysis of Stylolites and Faults
The coefficient of variation (C v ) is calculated as the ratio of the standard deviation of distances between features to the mean of distances between features (Gillespie et al., 1999) 10.1029/2019TC005965 where SD(s) represents the standard deviation of the sample of distances between each fault or stylolite and s represents the mean of the same sample distances. A C v of 0 represents a periodic distribution, whereas a C v of 1 represents a perfectly random distribution, and a C v greater than one represents clustered features (Cox & Lewis, 1966). Coefficient of variation analysis has been carried out on one-dimensional transects of discrete features to show the self-arrangement of stylolites (Railsback, 1998) and veins (Fagereng, 2011a;Gillespie et al., 1999). Following the approach of Gillespie et al. (1999), we threshold the distances between features at incrementally larger values to discern over what length-scale clustering is apparent. To compare stylolite and fault distributions within each lithological unit, we present these data as C v versus threshold thickness normalized to the thickness of the simplified lithological unit ( Figure 7).
Whole sample C v values for faults are smaller than for stylolites throughout the studied threshold range, representing less intense clustering of faults than stylolites ( Figure 7). This is also reflected in the unthresholded single value analysis (C stylolite v = 6.993, C fault v = 4.426). Stylolite spacing C v -threshold scaling relationships appear to follow a power-law relationship within Unit IV (Figure 7; cf. Gillespie et al., 1999). The degree of clustering of stylolites increases with smaller spacing in almost all lithologies; stylolites spaced >1 m are randomly to periodically distributed in all lithologies ( Figure 7). The periodic nature of stylolites spaced at thresholds >10 m probably reflects edge effects of the 350-m section of core studied, limiting analysis over this scale. This is also shown by low C v values in the shallowest calcareous mudstone, where only seven stylolites are present over a short section of core, compared to the entirety of Unit IV, where 873 stylolites are present (Figures 6 and 7). Stylolites are most clustered when spaced <1 cm apart in the deepest chalk unit. In all other considered lithologies, stylolites spaced 1 cm to 1 m are similarly clustered regardless of lithology ( Figure 7).
Fault spacing C v -threshold scaling relationships are more variable between lithologies than stylolites, though faults still cluster at spacings <1 m in randomly distributed horizons throughout Unit IV (Figures 6 and 7). Clustering within the marl, calcareous mudstone, and chalk lithologies is similar over all studied scales. Almost all faults spaced less than ∼1 m in these lithologies have C v > 1, indicating clustering. When faults are ∼1 m apart, C v is ∼1, indicating a random distribution (Figure 7). Faults in the debris flow are consistently more clustered than in other lithologies, maintaining a higher C v throughout almost the entire threshold range. In stark contrast, the 12 faults within the mudstone are periodic throughout the studied threshold range (Figure 7).

Tectonics 6. Strain Within Unit IV
Strain within Unit IV of Site U1520 was accommodated by mixed brittle and viscous deformation in the form of faults and stylolites. At sites of stylolite and fault co-occurrence, stylolites are mostly cut by fault surfaces, suggesting faulting occurred either more recently or more rapidly than stylolite growth. Faults and stylolites have therefore occurred over distinct time scales; stylolites appear to have accommodated strain more gradually through the history of the sediment while faulting occurred as a transient process, accommodating displacement rapidly before ceasing (Gratier et al., 1999). Unit IV can be traced landward and downdip to the active plate interface (Figure 1). With this in mind, we consider how vertical uniaxial strain from sedimentary load has been accommodated in Unit IV prior to subduction and describe how clustered faults and stylolites affect the initial rheology of sediments which could be accommodating strain at the plate interface.

Mass and Volume Loss Within Stylolites
Isocon plots use a known immobile species or, as in this case, an element present in a relatively insoluble mineral, to quantify mass changes between comparable analyses. Aluminum is present in clays and other relatively insoluble minerals (Aharonov & Katsman, 2009;Gratier et al., 2015;Hickman & Evans, 1995), commonly visible around and within calcareous fossils in the nearby protolith (Figures 3d-3g). Original depositional clay content of the analyzed areas within and without the stylolite is likely different as stylolites are more likely to form in areas of elevated clay concentration (Toussaint et al., 2018). The self-organizing nature of pressure solution, however, is likely to amplify these original differences as insoluble minerals are passively concentrated in zones of dissolution as more soluble carbonate is removed (Heald, 1955). We calculate bulk mass change by plotting the isocon through the plot origin and aluminum (Figure 8, Grant, 1986) and note that this method will overestimate mass loss if the stylolite formed in an area of initial high clay content. Titanium is also commonly used in isocon analyses, but low concentrations in our samples yield large errors and titanium was therefore not used here. From the bulk mass change, relative changes in calcium wt% are calculated (Figure 8).
Despite minor plagioclase (<10 wt%) throughout Unit IV (Figure 6), none was found in our EDS maps, reiterating that calcite is the dominant calcium-hosting phase losing mass on stylolites in originally calcite-rich sediments. The loss of calcite is consistent with (1) dissolution of fossils (Figures 3d-3g), (2) clusters of stylolites forming peaks of increased frequency with high CaCO 3 content (Figures 6 and 7), and (3) the relatively higher solubility of calcite, compared to silicates, under the studied conditions (Heald, 1955;Plummer & Busenberg, 1982), showing these data are consistent with mass loss by dissolution in calcareous rocks. We cannot constrain the amount of CO 2 lost from calcite in our samples. The bulk calcium mass change is therefore used to calculate volume change using calcite density and porosity.
Mass change in stylolite seams relative to the wall rock varies but is generally dominated by loss of Ca and gain of Mg, Si, Al, Ti, K, and Fe ( Figure 8). Bulk mass change calculated within stylolite seams shows loss of Ca-rich material, probably calcite, and concentration of aluminosilicates in highly localized zones (Figure 8). Relative total and calcium mass change estimates presented here are -11% to -45% and -3% to 49%, respectively ( Figure 8). Despite the possibility of overestimating mass loss due to initially higher concentrations of clay in stylolite-forming volumes, these mass losses are consistent with previous estimates of bulk mass change of -16% to -44% and CaCO 3 change of -15% to -45% from early pressure solution in marl from Lacroix et al. (2015).
Relative elemental mass changes in stylolites are variable between individual seams, but the majority of strain derives from porosity loss within the seam due to closer grain packing and fine aluminosilicates filling in intergranular voids (Figure 3). There is, again, the possibility that some of this calculated porosity change reflects initial spatial differences. We may overestimate the temporal volume change if stylolites formed on clay-rich, low-porosity horizons, or underestimate it if stylolites reflect initial high-porosity horizons where local normal stresses were elevated. We assume a porosity change of -40% based on porosity calculation from SEM BSE images (Figure 4), producing volume losses of 56% and 67-73% for bulk stylolites and individual seams, respectively ( Figure 8). Our calculated bulk volume losses from dissolution are similar to those estimated from deformed shales in the Shimanto Belt in Japan (13-54% Kawabata et al., 2007)

Strain in Stylolites
We now convert volume loss to shortening strain per individual stylolite and scale up our observations, using stylolite distribution, to determine the vertical shortening throughout Unit IV. As we are dealing exclusively with shortening strain, we express longitudinal strain values as positive for decreases in length throughout the remainder of this work to more clearly express the degree of shortening. Observational and chemical analyses of stylolites indicate uniaxial, vertical shortening by pressure solution, but quantification of this Numbers adjacent to elements are their multiplication factors. Dotted lines on isocon plots correspond to zero bulk mass change. Solid lines plotted through Al show actual bulk mass change, shown in percent in lower right of each plot. Horizontal error bars are 1 standard deviation of bulk protolith above and below stylolite; vertical error bars are 2σ from wt% analysis from EDS data. For details of mass change and strain calculations, see text. Plot (e) shows stylolite-hosted strain per recovered meter and stylolite frequency per recovered meter versus average CaCO 3 content within each core section. CaCO 3 content is from coulometric analysis performed during IODP Expedition 375. Filled points and lines correspond to the left axis and are colored according to depth downhole (see colorbar), colored points are strains per recovered meter calculated assuming a strain on each stylolite seam of 0.7 (e.g., a and b), whereas colored lines are strain ranges calculated assuming strains on stylolite seams from 0.1 to 0.9. Unfilled black circles are stylolite counts per recovered meter and correspond to the right axis. The exponential fit outline in part (e) corresponds to stylolite counts per recovered meter. As we are dealing solely with shortening strains, we show them as positive.

10.1029/2019TC005965
Tectonics strain is limited in scope to what is measurable within recovered core and by detailed study of selected, representative seams. We provide a limit for strain hosted by a singular stylolite seam by assuming all chemical volume loss is parallel to vertical stress (σ v ). We then use an idealized stylolite, composed of seven 50 μm wide seams of shortening 0.7 (a simplified approximation of the stylolite in Figure 8d), to calculate bulk uniaxial strain for recovered core within Unit IV. We present our data with error bars illustrating a range in calculated strain for each individual stylolite seam from 0.1 to 0.9 (Figure 8e). Bulk strain was calculated as change in length per recovered core length for each core section from holes U1520C and D using the equation: and ε recovered is the shortening strain in recovered core, L recovered is length of recovered core, L o is the undeformed thickness of recovered core, L styl is the thickness of stylolites within recovered core, calculated using the aforementioned idealized stylolite and the frequency of stylolites observed within recovered core, and ε styl is the shortening strain on an individual stylolite (0.7 from isocon plots and porosity loss; Figures 8a and 8b). Where core sections spanned the contact between two lithostratigraphic units, the percentage recovery was assumed to be the same either side of the boundary. The resultant shortening for each core section is plotted per recovered meter of core against average CaCO 3 content for the core section in Figure 8e. Also plotted is stylolite number per recovered meter within each core section, which shows an approximately exponential increase in strain as CaCO 3 contents increase above 70%. Increased stylolite-hosted strain with increasing CaCO 3 content is also particularly clear in the deeper sections of the core (Figure 8).

Controls on Fault and Stylolite Distribution
Both faults and stylolites cluster when spaced <1 m apart (Figure 7), meaning they are more likely to form near an existing feature than elsewhere (Gillespie et al., 1999). Clustering of features therefore suggests either: (1) feature formation changes conditions to make further nearby feature formation more likely, or (2) conditions in areas hosting many clustered features must be preferential for their formation. While some factors controlling the distribution of faults and stylolites are shared (differential stress, fluid pressure, etc), the distinct nature of faulting and pressure solution (i.e., physical vs. physiochemical) means that many of the factors controlling fault and stylolite distributions are distinct. Faults likely span several small-scale lithological boundaries with far greater lengths than stylolites and may have initiated in different lithologies before propagating through the core. We therefore refrain from suggesting any bulk controls on faulting based upon observations from drill core and focus on stylolites and stylolite-fault interaction.

Controls on Stylolite Formation
Because numerous factors influence stylolite formation (e.g., Toussaint et al., 2018), despite the large range of properties measured by IODP Expeditions 372 and 375 (Wallace et al., 2019), a peak in stylolite frequency does not always coincide with a measured change in core properties. We therefore do not attempt to fully characterize all factors which could influence stylolite formation but rather suggest some indicative variations which correlate with stylolite frequency.

Mineralogical Control on Dissolution
At a macroscopic scale, stylolites in Unit IV form anastomosing networks with low-amplitude, long wavelength roughness (Figure 3a). This is probably due to the lack of distinct mineralogical heterogeneities in a sediment dominated by carbonate and clay (Ebner et al., 2010;Railsback, 1993). Sediments in Unit IV rarely host more than 20% quartz and feldspar, and there are no discernible peaks in stylolite frequency correlating solely with an increase in quartz, feldspar, or clay content ( Figure 6). In contrast, stylolite frequency within clusters shows a correlation with CaCO 3 content throughout Unit IV (Figures 6 and 8e). This is logical as calcite, which the CaCO 3 probably represents, is more soluble than the other minerals present in Unit 10.1029/2019TC005965 Tectonics IV at the in situ P-T conditions (Heald, 1955;Plummer & Busenberg, 1982) and appears to be the main dissolved material (Figures 3d and 3e). Having said this, stylolites are not ubiquitous throughout all high CaCO 3 content lithologies, suggesting that high CaCO 3 content influences the intensity of dissolution rather than the location of stylolite formation. This CaCO 3 content control on dissolution magnitude likely exists between a lower bound where its presence allows dissolution and an upper bound where there is not enough phyllosilicates to form a stylolitic residue and localize pressure solution (cf. Zubtsov et al., 2004). This is reinforced by the total absence of stylolites in the mudstone, where CaCO 3 content is very low, the upper bound, as observed after dissolution is likely <95 wt% based upon observations of abundant stylolites within the CaCO 3 -dominated chalk at the base of Unit IV.

The Relationship of Stylolites and Lithology
Stylolites are most common in the deepest chalk unit, ∼650 out of 873 total stylolites occur there, and as CaCO 3 content is high throughout the chalks, this stylolite clustering requires an additional factor to explain heterogeneous stylolite development in space. On a submillimeter scale, stylolites dissolve relatively large (<200 μm) fossils throughout the sediments on discrete localized seams (Figures 3d and 3e). Fine clay leftover from dissolution is preserved within some stylolites as fine, elongate, grains with seam-parallel long axes (Figure 3e). Dissolution of a grain would continue until the grain contact area becomes so large that stress on that contact is below a critical value that is needed for pressure solution to operate, and grain boundary healing occurs (Rutter, 1983;van den Ende et al., 2019). Over the stress history of the sediments, grain contact surfaces probably dip below this critical value when the grain is elongated parallel to the seam, developing a shape-preferred orientation (Figures 3d and 3e). Such a shape-preferred orientation adjacent to stylolites has also been reported by many previous workers (Durney, 1972;Ebner et al., 2010;Toussaint et al., 2018). This interaction could explain the undulating nature of stylolites around and between fossils as they follow the path of maximum positive normal stress perturbation relative to the background (Figures 3d and  3e). Coupled with the accelerating effect of clay-rich and mineralogically diverse lithologies on dissolution (Aharonov & Katsman, 2009;Hickman & Evans, 1995), a clay-rich horizon with fossils interacting to concentrate stress would produce more rapid development of stylolites than a lithology without either clay or fossils. Indeed, Railsback (1993) shows that abundant contrast of clasts with a clay-rich matrix would host greater dissolution than either more abundant clasts or a more clay-poor matrix.
In contrast to fossiliferous lithologies, burrow margins within bioturbated lithologies host more well-developed stylolites and are more enriched in insoluble clay than both the internal volume of the burrow and the surrounding material ( Figure 3c). This enrichment suggests that where an initial concentration of clays or other insoluble material exists, dissolution is more localized and forms zones of well-developed stylolites, in agreement with observations from experimentally deformed diatomite and gypsum-illite composite by Gratier et al. (2015) and modeling by Aharonov and Katsman (2009) and Hickman and Evans (1995). Bioturbated and fossiliferous lithologies therefore have contrasting effects on stylolite development; bioturbated lithologies concentrate fewer, more intense, stylolites locally on burrow margins, while fossiliferous lithologies host many, less intense, more spaced stylolites.

Porosity Around Stylolites
Porosity is lower in stylolites than in the immediately adjacent material (Figure 4), suggesting porosity is lost by material dissolution causing closer grain packing (Macente et al., 2018;Pluymakers & Spiers, 2015). Larger porosity losses in more porous sediments (19R1W in Figure 4) show stylolites develop better in finer grained, higher porosity sediments. Similar to previous workers, we find zones of increased porosity immediately adjacent to stylolites (Figure 4, Heap et al., 2018). This is consistent with smaller grain contact areas supporting higher stresses, therefore localizing dissolution. Microfossils in stylolites are partially dissolved, showing porosity reduction where stylolites developed within sediments hosting fossil-bound void space (Figures 3 and 4). Absence of overgrowths adjacent to stylolites indicates that calcite did not precipitate immediately adjacent to the natural stylolites studied here, suggesting either precipitation rates were prohibitively low (i.e., precipitation-limited pressure solution), or there was limited solute available for precipitation. Diffusion-limited and precipitation-limited pressure solutions have been reported at low temperature (<150°C) for uniaxial compression of calcite at low (<0.04) and moderate strains (>0.04), respectively (Zhang et al., 2010). Although not well constrained, we raise the possibility of higher stylolite-marginal porosity promoting throughflow of fluids, allowing faster solute transport away from sites of dissolution and effectively limiting precipitation immediately adjacent to stylolites. Throughflow of fluids would likely 10.1029/2019TC005965 Tectonics freshen pore water at depth by lateral flow. This is not reflected in downhole chlorinity data at the specific depths where we observe stylolite clusters, but pore water sampling was too sparse to identify localized flow around individual stylolite clusters (Wallace et al., 2019). In concert with reduced grain size and slight shape-preferred orientation bordering stylolites (Figure 3), dissolution may also have occurred within a volume surrounding the stylolite, as previously suggested by Ebner et al. (2010). The outward propagation of the dissolution front from a clay-rich seam coincident with porosity reduction within stylolites has also been recognized by Macente et al. (2018) in experiments on a halite-biotite composite. This may explain how more well-developed stylolites grow to encompass several seams of phyllosilicates (Figure 3).

Temporal and Spatial Variations in a Stylolite Population
Given the lithological, physical, and chemical controls on the self-organising process of stylolite formation, the final distribution in Unit IV reflects spatial changes in carbonate and clay content, porosity and sedimentary structure, and grain-scale stress concentration and relative solubility (Figure 6), in decreasing order of importance. It seems the stylolite population within Unit IV is controlled dominantly by protolith properties rather than by stylolite formation altering initial conditions. The reasons for this are that stylolite seams (1) individually appear to have very localized regions of influence, only altering grain shape within tens of micrometers (Figures 3d and 3e), (2) form very closely without directly interacting, as would be expected if exerting a strong control on stylolite formation (Figures 3b and 3c), (3) correlate in frequency with peaks in bulk properties such as CaCO 3 content ( Figure 6), and (4) occur in frequency peaks which are randomly distributed over the length of Unit IV when spaced >1 m apart (Figure 7). It has been postulated that the concentration of clays during stylolite formation accelerates and localized further dissolution, meaning stylolite formation may cause further stylolite formation there (i.e., strain localization) (Aharonov & Katsman, 2009;Gratier et al., 2015;Hickman & Evans, 1995;Macente et al., 2018). This self-organization means that minor initial perturbations in the conditions outlined earlier are amplified once dissolution on stylolites begins. The subsequent temporal evolution of a stylolite is therefore conditional upon every previous feedback of evolving conditions with dissolution on that stylolite. It is entirely possible that when stylolites occur within tens of micrometers of one another, they are also influenced by changes in conditions resulting from adjacent stylolites, attested to by increasing clustering of closely spaced stylolites (Figure 7), but this is not a control on the macroscale distribution of stylolites.
Comparing stylolite frequency from 509 to 848 mbsf at Site U1520 to that from 790 to 1,350 mbsf on the Ontong Java Plateau to the north (Lind, 1993), it is clear a far higher frequency of less well-developed stylolites is present at shallower depths at Site U1520. Lind (1993) describes an erratic, highly localized, stylolite frequency, overall increasing in abundance with depth. Stylolites are said to be irregular in size and localize along burrows and clasts of contrasting mineralogy. This description is similar to what is seen at Site U1520, where randomly distributed peaks in stylolite frequency increase in intensity with depth ( Figure 6), burrows localize intense stylolite-rich regions along their margins (Figure 3c), and more fossiliferous regions host more spaced, less well-developed stylolites (Figure 3c). Carbonate content in the Ontong Java Plateau sediments is >90%, similar to the chalk at the base of Unit IV (Wallace et al., 2019). The difference in frequency is probably a result of a contrast in the intensity of stylolite development at each site; the limited depth range of Unit IV (509-848 mbsf) means that stylolites at Site U1520 are less well-developed than they would be at the depths studied by Lind (1993). Our counting is therefore more likely to be sensitive to less well developed stylolites, leading to a higher frequency. Regardless, controls on stylolite formation appear to be similar at both sites, suggesting that our observations of calcareous-pelagic sequences are valid beyond the Hikurangi sedimentary input sequence.

Interactions Between Stylolites and Faults
Microscale observations of faulting within Unit IV were limited, but faults appear to coincide with clay-rich weak layers prior to later alteration of their damage zone ( Figure 5). Thick (≤3 mm) clay-rich seams, concentrated by pressure solution, are present astride normal faults. These seams likely formed prior to fault formation, based on good preservation of the fault surface and the low strain rates at which pressure solution occurs relative to brittle fault formation by fracturing ( Figure 5, Gratier et al., 2013;Paterson & Wong, 2005). The brittle faulting therefore probably localized within pre-existing clay-rich seams, exploiting thin zones of relative weakness from previous viscous deformation (Viti et al., 2014). Both fault textural and distribution data, therefore, suggest stylolites are at least partially responsible for fault localization in Unit IV.

Tectonics
Where brittle fault surfaces are localized on strain-weakened low-angle surfaces, slip would result in centimeter-scale damage zones comprising both tensile and shear fractures adjacent to the slip surfaces ( Figure 5). These fractures represent new fluid pathways, rapidly increasing the surface area on which pressure solution can occur (Gratier et al., 1999). Strain rates on individual stylolite seams studied here (of strain 0.7; Figures 8a and 8b) range from 3.7 × 10 −16 to 2.2 × 10 −14 s −1 for ages of 60 and 1 Ma (Figure 2d), respectively. Although these are at the slow end of typical average geological strain rates in deforming zones (Fagereng & Biggs, 2019;Pfiffner & Ramsay, 1982). Our calculations of strain from stylolite frequency likely underestimate strain as it does not account for dissolution outside stylolite seams.

Modeling Uniaxial Strain Accommodation by Pressure Solution
The age, depth, and chemical constraints on the stylolite population described in this study provide a unique opportunity to compare strain accommodation from published constitutive models of pressure solution (Rutter, 1976;Spiers et al., 1990Spiers et al., , 2004Zhang et al., 2010) to strains calculated from direct observations of core from Site U1520 (section 6.2). The population of stylolites described here has been well characterized under well-defined P-T conditions, providing a directly comparable natural analog for the pressure solution described in the aforementioned models. We therefore use data from Site U1520 to construct a history of physical conditions and strain from intergranular pressure solution for the ∼60 Myr history of Unit IV.
As pressure solution is rate limited by the slowest of the serial processes of dissolution, diffusion, or precipitation, it is important to select the most appropriate model. The presence of partially dissolved fossils suggests dissolution is not the rate-limiting factor for pressure solution on stylolites (Figure 3), but as transitional diffusion or precipitation-limited behavior has been reported for calcite under conditions analogous to those studied here (Zhang et al., 2010), discerning the dominant rate-limiting factor from diffusion and precipitation is less clear. Based upon the observed lack of overgrowths adjacent to stylolites, we earlier inferred fluid flow and solute removal away from sites of dissolution in a system that is at least partially open (section 7.1.3), limiting local reprecipitation. If solute was rapidly removed from sites of dissolution, geochemical gradients would have been high at sites of dissolution and diffusion-limited pressure solution would have been faster than modeled in a closed system. We therefore use equations for diffusion-limited pressure solution as a lower bound on strain rate. As the number of assumptions within the modeling is large, we limit our discussion of absolute strain magnitudes, instead focusing on the location of pressure solution localization onto stylolites and which factors within the modeling could be better constrained.
We use the age model described earlier and reported in Wallace et al. (2019) to calculate past sediment depth at each depth and time step before calculating temperature and vertical stress using density data and a modern geotherm derived from measurements collected at the site ( Figure 2). We calculate stress conditions for pore fluid factors of 0.4 and 0.95 for comparison. Solubility (Plummer & Busenberg, 1982) and diffusivity (Nakashima, 1995) are calculated from estimated temperatures (20-33°C; section 3.2) and are used with a relatively thin grain boundary fluid film thickness of 1 × 10 −9 m (Renard et al., 1997) and an equation accounting for measured porosity at each depth (Pluymakers & Spiers, 2015) to calculate strain rate using published equations for diffusion-limited pressure solution (Rutter, 1976;Spiers et al., 1990Spiers et al., , 2004Zhang et al., 2010). We scale this modeled strain rate by the calcite content at that depth, both linearly and following the scaling relation of Zubtsov et al. (2004) (see section 8.2). Full details of the modeling methodology are available in Texts S1 and S2 in the supporting information.
Several assumptions are made within our history of physical conditions. Not least that data collected at a point during IODP Expedition 375 are representative for density, composition, porosity, and temperature gradient in the incoming sequence. We extrapolate these values back in time using an age model based upon macrofossils, microfossils, and palaeomagnetic reversal ages, all of which have varying degrees of dating error (Figure 2d). The only assumed constant condition throughout the history of the sediment is the modern temperature gradient from the seafloor downward, applied at each time step during the history of deposition. The main consequence of the temperature gradient is variability in the solubility of calcite, which is inversely correlated with temperature (Plummer & Busenberg, 1982). A higher temperature gradient, possibly due to higher heat flow from younger underlying oceanic crust (Villinger et al., 2002), would therefore reduce calcite solubility early in the sediment history. Therefore, our values probably overestimate 10.1029/2019TC005965 Tectonics calcite solubility in the early history of the sediment. Further to this, no attempt is made to quantify calcite solubility with increasing stress due to the lack of empirical published data.
Following constituitive equations for diffusion-limited pressure solution, grain-scale diffusivity effects within the model are captured within the diffusion coefficient and mass transfer width terms (Rutter, 1976(Rutter, , 1983. A simple Arrhenius temperature dependency (Nakashima, 1995) is unlikely to capture the complexity of electrochemical and grain boundary effects on the diffusion coefficient (D), while a single value for the mass transfer path width (S) ignores the complexity of subgrain size fluid pathways throughout grain boundaries between like and unlike minerals (Renard et al., 1997;Zubtsov et al., 2004). Furthermore, these values are likely to evolve with time and strain as well as the physical background conditions modeled here, possibly explaining the difference between our modeled and observed strains (Figure 11). We note that we have used a low-end value for the grain boundary fluid film thickness and current rather than initial porosities-these assumptions will tend to underestimate strain rates. Porosity within the sediments is present within and without spherical fossils ( Figure 4); meaning, porosity estimates from the core are likely greater than what is realistically available for precipitation of solute, if such precipitation does not occur within shells. The intrafossil porosity is not accounted for in the model, and the resulting strains are therefore likely to be overestimated. We consider this to be of relatively minor importance given the number of assumptions in the modeling and choose to focus our discussion on the distribution, rather than magnitude, of strain from pressure solution.

Modeling Results and Pressure Solution Controls
Modeled strain rates and cumulative strains are shown for the top and bottom of Unit IV in Figures 9e and 9f and for the whole of Unit IV in Figures 10c and 10d. Throughout the history of the sediment the dominant control on the modelled strain rate is the overlying sediment depth (Figure 10a). Sediment depth is controlled by the age model for the site (Figure 2d) and expressed in the modelling by effective normal stress (Figures 9d and 10a) and temperature (Figure 9b). After initial deposition, slow accumulation of sediment caused strain rate to increase gradually in the deepest 125 m of Unit IV during the first 50 Myr of the modeled history (from >10 −16.5 to ∼10 −16 s −1 ; Figures 2e, 9e, and 10c). These strain rates are relatively slow by geological standards (Fagereng & Biggs, 2019;Pfiffner & Ramsay, 1982) but were sustained long enough (until 15 Myr ago) to accommodate strains of up to ∼0.08 before the upper ∼215 m of Unit IV were deposited (Figures 9a, 9f, and 10e). Around 15 Myr ago, the age model dictates continuous sedimentation rate increases up to the present (Figure 9a). This manifests as fluctuating, relatively rapid increases in modeled strain rate (up to ∼10 −15.5 s −1 ) over ∼14 Myr, until rapid loading of turbidites in Units I-III in the last ∼1 Myr increases strain rates further (up to ∼10 −15 s −1 ; Figures 9e, 9f, 10c, and 10d). This strain rate increase has the most marked effect in the deepest 100 m of Unit IV, where sediments are over 25 Myr old (Figures 9f and 10d). At a given time-step modeled strain rates generally increase with depth, fluctuating locally because of local variations in porosity and calcite content (Figure 10d). Strain rates are highest between 825 and 850 mbsf, where calcite content is high, with local maxima coinciding with large fluctuations in measured porosity (e.g., ∼830 and ∼840 mbsf in Figure 10c). Unsurprisingly, this depth range of higher modeled strain rates, and the local maxima within it, has the highest modeled strain within Unit IV (Figure 10d). Final strains at the local maxima are up to 0.05 above the background within this depth range, even when the lower strain sediments are >10 Myr older (∼840 vs. ∼845 mbsf in Figure 10d).
The constant geothermal gradient applied during the modeling means temperature-dependent factors effectively become depth dependent, and the variation of solubility and diffusivity with temperature can be discussed alongside that of depth-dependent factors such as effective normal stress. During the first 50 Myr of our model, the deepest calcareous sediments (>700 mbsf) maintain low effective normal stresses (<2 MPa) and temperatures (<10°C), corresponding to high solubilities (>10 −5.6 m 3 m −3 ; Figures 9b-9d, 10a, and  10b). The strain accumulated over this period (up to 0.08) is not negligible compared to final totals in our model (up to 0.15), but rapid loading in the more recent history of the model increases strain rates from ∼10 −16 s −1 to around 10 −15 s −1 , adding an extra strain of 0.07 over the final 15 Myr of Unit IV at average strain rates of approximately 10 −15.6 s −1 . This is similar to the accumulated strain over the previous ∼50 Myr and, while not surprising, illustrates how a nonlinear history of sedimentation or stress loading may affect rates of pressure solution during the history of deformation. Changes in effective normal stress, reflected in the modeled strain rate (Figure 10c), are even more apparent when pore fluid pressure is 10.1029/2019TC005965 Tectonics increased to near-lithostatic levels, reducing the stress driving dissolution and the modeled strain rate, and limiting the resultant strain ( Figure S1).

Comparing Modeled and Calculated Strains from Pressure Solution in Unit IV
Our calculation of stylolite-hosted strain within Unit IV is based on stylolite frequency and strain on each individual stylolite seam (Figure 8), meaning our total calculated strain reflects stylolite frequency in core image logs (Figures 6 and 10). This stylolite-hosted strain is likely a minimum value for strain within the whole of Unit IV, as pressure solution is likely to act upon much of the depositional layering as a diagenetic process (Gratier et al., 2015), and we have not quantified strains caused by granular flow or faulting. Modeled strains, however, are calculated using equations derived from experiments of distributed dissolution and bulk shortening rather than that localized on stylolites (Rutter, 1976(Rutter, , 1983Zhang et al., 2002Zhang et al., , 2010. Strains are therefore comparable within the depth range of stylolite clusters or over length scales encompassing several clusters (>1 m). Modeled strain from pressure solution at hydrostatic fluid pressure is almost always higher than observed strain, each with peaks which rarely coincide (Figure 11a). When studied at the length scales of lithological units (10s to 100s m) or core sections (∼10 m), calculated and modeled strains rarely agree, with modeled strains almost always exceeding the measured magnitude of strain (Figures 11b and 11c). At near-lithostatic pore fluid pressures, the model reproduces the average strain magnitude when compared over larger length scales (Figures 11e and 11f) but does not reproduce observed stylolite clustering, giving lower than measured strain within stylolite clusters and higher than measured strain between clusters (Figure 11d). As modeled strain uses bulk properties throughout the core and stylolites host highly localized dissolution, it seems that stylolite formation is controlled by local perturbations not accounted for in our models. This is not surprising as we cannot model unknown initial heterogeneities  Plummer and Busenberg (1982), (c) lithostatic and hydrostatic stress, (d) effective normal stress, (e) strain rate from diffusion controlled dissolution of calcite according to published models for pressure solution (Rutter, 1976;1983),

10.1029/2019TC005965
Tectonics in the sediment column, and our models are therefore averages over multiple clusters. Despite our attempt to integrate the evolution of some of the parameters controlling stylolite evolution, the use of the constitutive equations of pressure solution (Gratier et al., 2013;Rutter, 1976Rutter, , 1983Spiers et al., 1990Spiers et al., , 2004 conceived and calibrated with constant or controlled parameter values does not accurately predict the location or magnitude of cumulative strain of intensive dissolution on stylolites. Beside the possible uncertainty on the values of the various parameters of the equations, the evolution of the parameter values with stylolite evolution is not yet enough constrained within our model, particularly the effects of self-organized processes. Porosity is variable both within Unit IV and locally around stylolites (section 7.1.3 and Figure 4). The coincidence of observed stylolite clusters and localized porosity peaks (e.g., at ∼760 and ∼805 mbsf) is consistent with observations of greater volume loss and increased resultant shortening strain in higher porosity zones (19R1W in Figure 4). Modeled strains have local high-frequency variations with depth due to variable porosity and calcite content measured throughout Unit IV (Figures 6c and 6d). With depth, some peaks in cumulative modeled strain coincide with peaks in stylolite-hosted strain (Figures 11a and 11d). Previous studies that show strain increases with increased contact stress over reduced grain contact area (Niemeijer et al., 2002;Pluymakers & Spiers, 2015;van-den-Ende et al., 2018;Zhang et al., 2010), suggesting that applied variations in strain rate with porosity may be, at least partially, representative of variations in pressure solution within natural basins containing calcareous sediment.
Composition is another important difference between the sediments studied here and the materials used to study the dissolution of pure calcite (Rutter, 1976(Rutter, , 1983Zhang et al., 2002Zhang et al., , 2010. Our attempt to account for the heterogeneous mineralogy of Unit IV by limiting strain from dissolution to soluble grains does not account for (1) the accelerating effects of electrochemical potential on pressure solution across mineralogically heterogeneous grain boundaries or (2) stylolite localization associated with increased clay fraction at particular horizons (Aharonov & Katsman, 2009;Greene et al., 2009;Hickman & Evans, 1995). Furthermore, stylolites are typical of pressure solution in dominantly calcareous sediments, whereas

10.1029/2019TC005965
Tectonics pressure solution within more clay-rich sediments would likely occur on depositional layering, accommodating vertical shortening without forming stylolites (Gratier et al., 2015). Evidence for this pressure solution in more clay-rich lithologies would be difficult to observe and has not been counted here as a straining feature; meaning, our measurements of stylolite-hosted strain are likely a minimum value for the whole of Unit IV. Previous experiments on pressure solution in two-phase mixtures found a linear scaling relationship exists between strain rate and soluble species content up to 50%, where it plateaus until 70% soluble species, before decreasing linearly at soluble species contents of 70-100% (Zubtsov et al., 2004). This was thought to be because diffusion in a fluid phase trapped along grain contacts between different minerals is faster than along grain contacts between the same mineral. The mineral assemblage varies significantly throughout Unit IV ( Figure 6) and variations in grain-scale factors such as these could account for variations in observed strain that have not been recreated in our model. We apply the scaling relation of Zubtsov et al. (2004) but find it does not better recreate calculated stylolite-hosted strain distribution or magnitude ( Figure S2). Better quantitative constraints to measure and model the effect of grain-scale processes on macroscopic strain, by both distributed and localized pressure solution, would aid the ability of models to recreate observed trends in strain.

Tectonics
Cumulative strains modeled over the history of Unit IV vary significantly in magnitude from those estimated from stylolite mass loss and frequency. Model results for grain sizes of 10 and 100 μm vary over three to four orders of magnitude, and agreement of calculated and modeled strains is only achieved at hydrostatic fluid pressures when the modeled grain size is ≥100 μm (Figures 11b and 11c) and at near-lithostatic fluid pressures when the modeled grain size is ≥50 μm (Figures 11e and 11f). The mean grain size measured in the chalks in Unit IV is <40 μm, with larger microfossils ≤200 μm (Figure 3), suggesting that either (1) the model overestimates strain rate from grain dissolution or (2) near-lithostatic pore fluid pressures have been sustained for long periods of the sediment history and (3) not all pressure solution forms stylolites or (4) strain within Unit IV is controlled by the dissolution of grains larger than those comprising the matrix (Figures 10  and 11). The dissolution of microfossils seen in microstructural observations suggests that the latter may be the case, though matrix grains also appear to be dissolved, showing a shape-preferred orientation adjacent to stylolites (Figure 3). Alternatively, the fine-grained nature of the sediment, and possible fluid compartmentalization by overlying low permeability mudstones (Figure 1b), may locally lead to greater than hydrostatic pore fluid pressures during the history of the sediment, reducing effective vertical stress and subsequent dissolution. However, there is no other evidence, such as fractures, veins, or fluid escape structures, indicative of sustained near-lithostatic fluid overpressures at Site U1520 (Wallace et al., 2019). A further point for consideration is that the diffusive distance (d) is assumed to be equal to the grain size in our model. It may be that in systems where lateral long-range solute transport may be operating (e.g., in an open system with advective transport of calcite) diffusive distance can be larger than grain size (Dewers & Ortoleva, 1990;Gundersen et al., 2002;Lehner, 1995).

Forward Modeling Pressure Solution to the Plate Interface
Having modeled pressure solution in Unit IV through the last ∼60 Myr to the present day, we also extend the dissolution model to show the final effect of pressure solution when the sediments reach the toe of the Hikurangi plate interface. We use the thickness of sediment overlying Unit IV from Site U1520 to the toe of the thrust on seismic line 05CM-04 to calculate vertical stress into the future (Figure 1 Barnes et al., 2020). Strain rate from pressure solution is calculated as previously.
The thickness of sediment added is less than 500 m at the toe (Figure 1), reflected by a change in normal stress of ∼4 MPa ( Figure S3). Similarly, a minor change in calcite solubility results from higher temperature with increased depth. The distance between Site U1520 and the toe of the megathrust is relatively short (∼20 km), occurring rapidly compared to the previous model and hosting shortening strain increases of less than 0.025 throughout Unit IV ( Figure S3). The small difference in strains indicates that the majority of pre-subduction dissolution-hosted strain has already been accommodated within Unit IV over the last (∼60 Myr) so that stylolites observed at Site U1520 are likely very similar to those existing within Unit IV at the toe of the megathrust.

Future Subduction of Unit IV and Rheology of the Current Plate Interface
Modeled and observed strains from pressure solution in Unit IV at site U1520 are comparable when averaged over tens of meters and using a diffusion length scale comparable to the coarsest sediment fraction. Forward modeling to the toe of the megathrust shows only a slight change in strain, suggesting the majority of pre-subduction compaction had already been accommodated by the time the sediment reached Site U1520 ( Figure S3). Depth-converted seismic profiles indicate the lower part of Unit IV is likely the protolith for the localized plate interface, where it is likely to accommodate shear strain during subduction (Barnes et al., 2020).

Rheological Variation Throughout Unit IV
Unit IV is highly variable with regard to mineralogy, porosity, and lithological texture ( Figure 6). These factors, along with partial lithification and abundant pore water, are likely to influence slip and deformation style and contribute to a heterogeneous plate interface (Barnes et al., 2020;Delle Piane et al., 2017;Mittempergher et al., 2018;Sample, 1990;Skarbek & Rempel, 2017;Wallace et al., 2019). Compounding this, the calcareous-pelagics in Unit IV overlie a sequence of volcanic conglomerates, completely rheologically distinct and even more locally variable (Barnes et al., 2020;Wallace et al., 2019). The rheology of intact lithologies comprising Unit IV at the plate interface would therefore already be locally variable, but pre-subduction 10.1029/2019TC005965 Tectonics strain on faults and stylolites will cause local changes in rheology from those predicted for the intact sediments of Unit IV.
Strains calculated from stylolites within Unit IV are not large. Uniaxial shortening strain for the entirety of Unit IV is ∼0.004 but represents local accumulation of insoluble clays within each discrete <1-mm-thick stylolite (Figure 3). Each stylolite does not host enough insoluble clay to significantly alter the rheology of the bulk rock, but where stylolites cluster the accumulation of clays within stylolites may become locally significant ( Figure 6). Stylolites cluster when less than 1 m apart (Figure 7), with clusters hosting up to 80 stylolites per meter (Figure 6). Clay concentrated on stylolites in these clusters represents strain-weakened horizons from early viscous deformation due to the weakness of clay relative to the encompassing carbonate (Boulton et al., 2019;Rabinowitz et al., 2018), possibly significant enough to locally alter bulk deformation style. Indeed, some faults in Unit IV are localized along clay-rich seams of similar thickness to those within stylolites (Figures 3 and 5). Clusters of stylolites increase in intensity with depth, increasingly altering lithological texture and resultant deformation style on increasingly interconnected anastomosing stylolite seams (Figures , 3, 6, and 12).
Faults within Unit IV are rheologically variable, with weak surfaces between two moderately weak zones of damaged material (Figures 5 and 12). Individual fault and fracture zones are more commonly considered than stylolites as zones of strain weakening, likely to host later deformation (Ferraro et al., 2018;Leah et al., 2018;Mitchell & Faulkner, 2009), but clustering of faults in Unit IV makes contemporaneous reactivation of several faults within a centimeter-to meter-thick volume likely at the plate interface (Figures 6, 7, and 12). In addition, pressure solution has concentrated phyllosilicates on fractures in damage zones, concentrating clays and further weakening areas hosting fault clusters (Figures 5 and 12).

Future Décollement Location and Deformation Localization
Within Unit IV, clusters of clay-rich stylolite seams or fracture-bound faults represent these low-volume fraction phases with lower shear strengths than the host rock, likely to localize deformation at the plate interface. Rheological heterogeneity is another key factor in deformation localization during reactivation of pre-existing weaknesses (Imber et al., 2008;Leah et al., 2018;Lyakhovsky et al., 1997;Willemse et al., 1997), and stylolite and fault clusters represent heterogeneous intervals throughout Unit IV. One particular interval, near the base of Unit IV (835-845 mbsf; Figure 6), hosts nearly 400 stylolites in two large clusters, increasing the volume fraction of clay present by localized pressure solution ( Figure 6). This increased clay fraction creates a shear strength contrast within the sediments, making strain more likely to localize there when subject to shear stress. Shallower in the sediment stack, mudstones within Unit IV (e.g., ∼720-740 mbsf) may also be significantly weaker layers likely to localize shear strain. The localization of deformation likely to occur upon reaching the plate interface would manifest as incrementally greater shear strain is accommodated in the weaker horizons rather than in the stronger intervals, creating a thin cumulative volume hosting high shear strain, including a likely horizon near the base of Unit IV ( Figure 12).

Deformation Mechanisms at the Plate Interface: Implications for Slow Slip
Slip on faults and fractures within their damage zones will be frictional, with properties controlled by that of the material on the fault. Where faults have localized upon pre-existing stylolites, therefore, slip character is likely to be controlled by phyllosilicates localized in the stylolite seam ( Figure 5). Where faults have formed away from stylolites, frictional slip is likely to reflect the bulk protolith, dominantly calcite ( Figure 6). Calcite and clays are likely to mix during slip, possibly by weak phase smearing (Rutter et al., 2013). Friction on these surfaces would likely be dominated by phyllosilicate-hosted slip with a lesser component of calcite. This would lead to a mixture of velocity strengthening and weakening conditions over varying conditions and could provide appropriate conditions for unstable slow slip on the plate interface (Boulton et al., 2019;Rabinowitz et al., 2018;Tesei et al., 2014).
Localized slip surfaces are also likely to form on weak, clay-rich stylolites, where they could accommodate shear strain by a combined frictional-viscous flow mechanism Gratier et al., 2013;Willemse et al., 1997). This entails frictional slip on an anastomosing phyllosilicate foliation with removal of material from grains within the foliation by pressure solution, similar to deformation modeled for quartz-bearing gouge by den Hartog and  or less phyllosilicate-rich halite or quartz-bearing gouge by Bos and Spiers (2002). As calcite solubility is high at low temperatures (Plummer & Busenberg, 10.1029/2019TC005965 Tectonics 1982), abundant calcite within these deforming zones will likely deform by pressure solution at the shallow plate interface. Intensive pressure solution would cause volume loss from the volume accommodating shear between clasts, possibly increasing relative clay volume fraction within the sliding portion of the plate interface shear zone and altering the shape of clasts (Gratier et al., 2015), thereby causing strain weakening. Alongside this, carbonates experience significantly more variable deformation at relatively low temperatures than their siliciclastic counterparts, including recrystallization (Kennedy & White, 2001), crystal plasticity (Verberne et al., 2013), and variable frictional properties (Verberne et al., 2014).
Within a plate interface shear zone hosting grain size reduction of heterogeneous materials over varying depths, pore fluid pressures, and temperatures (Barnes et al., 2020), it is likely that pressure solution will be locally variable on the submillimeter scale, removing and redepositing material locally around grains or clasts. At shallow depths, the temperature gradient of the northern Hikurangi Margin is thought to be relatively high (Antriasian et al., 2018), corresponding to a rapid reduction in calcite solubility as materials are buried (Figure 12). The ability of pressure solution to cause removal of calcite, therefore, is greatly reduced outside regions of low temperature, likely corresponding to less than ∼150°C or 10-to 15-km depth ( Figure 12 Antriasian et al., 2018). This is a similar depth to the deeper end of the slow-slipping zone, possibly indicating that pressure solution of calcite may decrease in relative importance as a mechanism accommodating deformation at the plate interface outside the area recognized to host slow slip.

Conclusions
Unit IV at Site U1520 hosts clustered faults and stylolites within muddy pelagic carbonates with variable composition. The frequency of stylolites increases with depth, whereas faults are more common within several shallower horizons. Stylolite frequency increases approximately exponentially with CaCO 3 content. Figure 12. Schematic overview of strain localization and weakening in Unit IV at Site U1520 and its possible future rheological behavior during subduction. Sketches with notes (above) show locally important textures and rheological notes from Unit IV at various stages during subduction (middle). Thermal gradient (below) is from Antriasian et al. (2018); resultant calcite solubility is calculated using the method of Plummer and Busenberg (1982).

10.1029/2019TC005965
Tectonics Stylolites comprise anastomosing individual seams, each hosting mass losses of 11-45% within a horizon of lower porosity compared to the host rock. Uniaxial shortening, accounting for porosity change, is high on individual stylolite seams (∼0.5) but small over greater distances because localized horizons host dissolution. Stylolites concentrate clays by dissolution, forming localized clay-rich surfaces within a mostly intact host rock. Faults in Unit IV are present in more lithified units, locally exploiting pre-existing stylolites. Slip surfaces host adjacent centimeter-scale damage zones which have themselves undergone pressure solution, showing that mixed brittle and viscous deformation occurs at very shallow depths in calcareous sediments. Stylolite distribution reflects an interplay of grain-scale physical and larger scale compositional factors, forming most commonly where calcite content is high, clay is present to accelerate dissolution, or high initial porosities localize dissolution.
We model pressure solution using published models for intergranular pressure solution of carbonates over the depositional history of Unit IV. To do this, we determine an absolute P-T-t history using porosity, density, age, and temperature constraints from IODP Expedition 375. Despite incorporating modern day composition and the temporal evolution of stress and temperature into the pressure solution model, the locations and magnitudes of strain within stylolite clusters observed within Unit IV are not well reproduced at hydrostatic or near-lithostatic fluid pressures. At hydrostatic fluid pressures, we find agreement at grain sizes corresponding to coarse fossils within the sediment, possibly consistent with porosity loss inferred from stylolite microstructure. At near-lithostatic fluid pressure, observed and modeled strains agree at more common grain sizes in the sample, though this requires sustained periods of high pore fluid during the ∼60-Myr history of Unit IV. Strain from localized dissolution on stylolites is likely overestimated using a bulk pressure solution model due to the complexity and evolution of subgrain-scale fluid pathways and diffusivities not fully captured within the applied model. We show that a linear scaling of calcite volume percent and strain rate is overly simplistic and that the relationship of grain-scale effects within pressure solution models must be better constrained with composition, time, and strain. We extend the model to the toe of the megathrust and show that dissolution causes negligible change in strain between present day and initial subduction of the sediment column at Site U1520.
Shear within the carbonaceous Unit IV at the plate interface is likely to localize in weak horizontal horizons including clay-rich stylolite clusters at the bottom of Unit IV, slipping by frictional-viscous mechanisms with variable frictional stability due to composition heterogeneity within stylolite clusters. Pressure solution of calcite is likely an important viscous mechanism acting on the grain scale throughout the slow-slipping zone and at shallow depths on the plate interface. Indeed, the down-dip end of slow slip events (15-20 km) roughly correlates with where carbonate solubility becomes negligible. Throughout the slow-slipping zone and the shallow plate interface, other deformation mechanisms in the chalks and clay-rich muds are likely highly variable. The role of carbonate in accommodating shear at the Hikurangi Margin is emphasized, as this dramatically alters slip characteristics and deformation style compared to siliciclastic sediment.