Burial-Related Compaction Modifies Intrusion-Induced Forced Folds: Implications for Reconciling Roof Uplift Mechanisms Using Seismic Reflection Data

Space for shallow-level sills and laccoliths is commonly generated by bending and uplift of overlying rock and sediment. This so-called ‘roof uplift’ produces forced folds, the shape and amplitude of which reflect the geometry of underlying intrusions. The surface expression of forced folds can therefore be inverted to constrain intruding magma body properties, whilst ancient forced folds provide a record of sill and laccolith emplacement. Deciphering how shallow-level intrusion translates into roof uplift is thus critical to enhancing our understanding and forecasting of magma emplacement. To-date, emplacement models and surface deformation inversions are underpinned by the consideration that roof uplift is, to a first-order, an elastic process. However, several studies have suggested inelastic processes can accommodate significant magma volumes, implying first-order roof uplift may be a function of elastic and inelastic deformation. In particular, seismic reflection images of forced folds above ancient sills and laccoliths have been used to argue that final fold amplitudes can be substantially less (by up to 85%) than the underlying intrusion thickness. Although these seismic-based observations imply elastic and inelastic deformation accommodated intrusion, these studies do not consider whether burial-related compaction has reduced the original fold amplitude. Here, we use geological (e.g. lithology) and geophysical (e.g. seismic velocity) information from the Resolution-1 borehole offshore eastern New Zealand, which intersects a forced fold and upper ~50 m of a sill imaged in 2D seismic reflection data, to decompact the folded sequence and recover its original geometry. We show the Resolution Sill is likely ~117–187 m thick, depending on the interval velocity for the entire intrusion, whereas the forced fold has an apparent maximum amplitude of ~127 m, corresponding to a sill thickness-fold amplitude discrepancy of up to 47%. Decompaction indicates the original maximum forced fold amplitude likely ranged from ~131–185 m, suggesting post-emplacement, burial-related compaction of this and other forced folds may be the source of apparent discrepancies between fold amplitude and intrusion thickness. Whilst seismic reflection data can provide fundamental insights into how shallow-level emplacement translates into roof uplift and ground displacement, we show decompaction and backstripping are required to recover the original fold geometry.

Space for shallow-level sills and laccoliths is commonly generated by bending and uplift of overlying rock and sediment. This so-called "roof uplift" produces forced folds, the shape and amplitude of which reflect the geometry of underlying intrusions. The surface expression of forced folds can therefore be inverted to constrain intruding magma body properties, whilst ancient forced folds provide a record of sill and laccolith emplacement. Deciphering how shallow-level intrusion translates into roof uplift is thus critical to enhancing our understanding and forecasting of magma emplacement. To-date, emplacement models and surface deformation inversions are underpinned by the consideration that roof uplift is, to a first-order, an elastic process. However, several studies have suggested inelastic processes can accommodate significant magma volumes, implying first-order roof uplift may be a function of elastic and inelastic deformation. In particular, seismic reflection images of forced folds above ancient sills and laccoliths have been used to argue that final fold amplitudes can be substantially less (by up to 85%) than the underlying intrusion thickness. Although these seismic-based observations imply elastic and inelastic deformation accommodated intrusion, these studies do not consider whether burial-related compaction has reduced the original fold amplitude. Here, we use geological (e.g., lithology) and geophysical (e.g., seismic velocity) information from the Resolution-1 borehole offshore eastern New Zealand, which intersects a forced fold and upper ∼50 m of a sill imaged in 2D seismic reflection data, to decompact the folded sequence and recover its original geometry. We show the Resolution Sill is likely ∼117-187 m thick, depending on the interval velocity for the entire intrusion, whereas the forced fold has an apparent maximum amplitude of ∼127 m, corresponding to a sill thickness-fold amplitude discrepancy of up to 47%. Decompaction indicates the original maximum forced fold amplitude likely ranged from ∼131-185 m, suggesting post-emplacement, burial-related compaction of this and other forced folds may be the source of apparent discrepancies between fold amplitude and intrusion thickness. Whilst seismic reflection data can provide fundamental insights into how shallow-level emplacement translates into roof uplift and ground displacement, we show decompaction and backstripping are required to recover the original fold geometry. To assess the relative importance of elastic and inelastic space-making processes during the formation of seismically imaged sills and forced folds, we demonstrate that our method should be applied to remove any post-emplacement, burial-related compaction signature.
Seismic reflection data reveal the current maximum amplitude (F max ) of buried forced folds can be up to 85% less than the measured maximum thickness (T max ) of underlying, crystallized sills or laccoliths (i.e., F max /T max < 1; Figure 1A) (Hansen and Cartwright, 2006;Jackson et al., 2013;Magee et al., 2013a). Such discrepancies between fold amplitude and intrusion thickness, particularly where F max /T max << 1, have been suggested to relate to the accommodation of magma by both elastic and inelastic deformation (Jackson et al., 2013;Magee et al., 2013aMagee et al., , 2017. Syn-intrusion, fracture-driven porosity reduction, faulting, and fluidization of the host rock around exposed sills confirms that inelastic deformation can partly and, perhaps in some instances, fully accommodate magma emplacement (Figures 1B,C) (e.g., Johnson and Pollard, 1973;Morgan et al., 2008;Schofield et al., 2012Schofield et al., , 2014Jackson et al., 2013;Spacapan et al., 2016). It has also been suggested that inelastic ductile strain and vertical compaction of deforming strata can cause fold amplitudes to decay upwards, particularly if D/d is <4 (Hansen and Cartwright, 2006;Jackson et al., 2013). Seismic and field data therefore provide evidence for the accommodation of magma by elastic and inelastic deformation, challenging the assumption that emplacement models need only account for elastic processes (e.g., Galland and Scheibert, 2013;Magee et al., 2013a;Holohan et al., 2017;Scheibert et al., 2017;Gerbault et al., 2018).
Seismic reflection data capture the current, and not necessarily the original, geometry of ancient intrusions and forced folds. For example, original fold amplitudes and sill thicknesses, and the ratio between them, may by modified by the: (i) migration of magma away from the seismically resolved intrusion (T max < T0 max ) coupled with little or no fold subsidence (F max > T max ) (e.g., Reeves et al., 2018); (ii) deflation of the sill in response to crystallization of and/or volatile release from the magma (T max < T0 max ; e.g., Caricchi et al., 2014), which could promote disproportionate fold subsidence (F max > T max ); (iii) erosion of the fold crest (F max < T max ) (Hansen and Cartwright, 2006;Jackson et al., 2013); (iv) interference with neighboring folds (Jackson et al., 2013), which could locally inhibit or enhance folding; and/or (v) burial and compaction of the folded sequence (F max < T max ) (Jackson et al., 2013). No study has yet quantified how post-emplacement, burialrelated compaction can modify forced fold geometries and amplitudes. Without incorporating an assessment of how burialrelated compaction has affected the seismically resolved forced fold geometry and amplitude, the role of inelastic processes in accommodating magma cannot be determined from seismic reflection data alone.
Here, we examine a saucer-shaped sill, the Resolution Sill, and overlying forced fold imaged in 2D seismic reflection data from the Canterbury Basin, offshore eastern New Zealand and intersected by the Resolution-1 borehole (Figure 2). The borehole penetrates the upper ∼50 m of the saucer-shaped sill, which can broadly be categorized as a gabbro. Using velocity information from Resolution-1 we aim to depth convert and decompact the seismic reflection data, allowing us to constrain the original maximum fold amplitude (i.e., F0 max ) and assess how burial-related compaction impacts fold geometry. We show that burial-related compaction modifies ancient intrusion-induced forced folds within sedimentary basins, reducing discrepancies between fold amplitude and sill thickness. Before using seismicbased examples of ancient intrusion and forced fold pairs to FIGURE 1 | (A) Seismically imaged maximum forced fold amplitudes (F max ) plotted against maximum measured thicknesses (T max ) of underlying sills or laccoliths from data within the: (i) Bight Basin, offshore southern Australia (Jackson et al., 2013); (ii) Exmouth Sub-basin, offshore north-western Australia (Magee et al., 2013a); and (iii) Rockall Basin, NE Atlantic (Hansen and Cartwright, 2006). See Jackson et al. (2013) and Hansen and Cartwright (2006) for information on error bars. (B) Field photograph showing folding of sandstone beds above the Trachyte Mesa intrusion in the Henry Mountains, Utah, USA. (C) Sketch showing changes in thickness of a massive red sandstone bed, shown in (B), over the Trachyte Mesa intrusion, which corresponds to a reduction in porosity (after Morgan et al., 2008). The best-fit curve is third order polynomial with an R 2 value of 0.49 (Morgan et al., 2008). Frontiers in Earth Science | www.frontiersin.org postulate emplacement mechanics at active volcanoes, it is essential to first account for burial-related compaction.

GEOLOGICAL SETTING
The Canterbury Basin spans onshore and offshore SE New Zealand and formed during Late Albian-to-Early Campanian rifting between New Zealand, Antarctica, and Australia (Figure 2) (Fulthorpe et al., 1996;Lu and Fulthorpe, 2004). The basement typically corresponds to greywacke and argillite metasedimentary rocks of the Torlesse Supergroup (Permian-to-Early Cretaceous; Figure 3) (Uruski, 2010). In the north of the basin, within the study area, syn-rift sedimentary strata deposited within graben and half-graben are dominated by the paralic coal measures of the Broken River Formation, and marine siltstones and mudstones of the Conway Formation (Figure 3) (Carter, 1988;Killops et al., 1997;Schiøler et al., 2011). Onset of post-rift, thermal subsidence in the Maastrichtian led to the deposition of the high-energy marine Charteris Bay Sandstone (Lower Paleocene), which is overlain by tuffs of the View Hill Volcanics, mudstones of the Conway Formation, and calcareous marine mudstones of the Ashley Formation (Figure 3) (Carter, 1988;Killops et al., 1997;Schiøler et al., 2011). Micritic limestones attributed to the Amuri Formation were deposited between the Early Oligocene and Early Miocene, although the majority of this time period corresponded to the development of a regional unconformity across much of the Canterbury Basin (Figure 3) (Carter, 1988;Killops et al., 1997;Schiøler et al., 2011). Uplift along the Alpine Fault, and an associated increase in the supply of terrigenous silt and sand, resulted in the deposition of the marine Tokama Siltstone, which locally contains tuffs belonging to the Harper Hills Basalt (K-Ar ages of 13.5 ± 0.4-11.0 ± 0.3 Ma), and overlying Kowai Formation (Early Miocene-to-Recent; Figure 3) (Sewell and Gibson, 1988;Lu et al., 2005).
Several discrete phases of intra-plate, post-Cretaceous magmatism and volcanism have been recorded in the Canterbury Basin, including the View Hill Volcanics and the Harper Hills Basalt (Figure 3) (e.g., Timm et al., 2010;Reeves et al., 2018). It has been suggested that volcanism occurred in response to decompression melting of upwelling heterogeneous asthenospheric mantle following localized removal of gravitationally unstable lithospheric material (Timm et al., 2010).

DATASET Borehole Data
Resolution-1 is located ∼50 km south of Christchurch (Figure 2) and was drilled in 1975 for Shell BP Todd Canterbury Services Ltd (Milne et al., 1975). The borehole was drilled in a water depth of 64 m and extends to a total depth of 1,963.05 m, intersecting the Resolution Sill between 1,911.5 and 1,963.05 m (Milne et al., 1975). Data available for the borehole include (Milne et al., 1975): (i) a well-completion report containing petrological descriptions of cuttings and sidewall core, sampled every 5 m between 1,910 and 1,958 m, and continuous core collected between 1,958.2 and 1,963.05 m within the sill; (ii) sonic ( T), gamma ray (GR), caliper (CAL), and spontaneous potential (SP) logs ( Figure 4); (iii) a petrophysical summary log plot; (iv) well-formation tops, ages, and lithological descriptions; and (v) K-Ar ages of 12 ± 2 Ma for the sill. Density logs, neutron porosity logs, thin sections, or photomicrographs are not available to corroborate the petrographic descriptions.
Resolution-1 has sparse time-depth information. To facilitate depth-conversion of the seismic reflection data, we therefore derived a time-depth curve by integrating sonic log data after using a median filter with a window of five samples to remove spikes caused by sample skipping (Figure 4). The sonic log data were also used to calculate a compressional wave (V p ) velocity log by taking the reciprocal of the interval transit time log and converting from feet to meters, and to define average interval velocities for different units (Figure 4). For example, the average interval velocity within the Resolution Sill intersected by the borehole is 5.2 km s −1 (Figure 4), with a standard deviation of 0.3 km s −1 . Although the average interval velocity of the sill where it is intersected by Resolution-1 can be defined (i.e., 5.2 km s −1 ), the borehole does not extend through the entire intrusive body; as a result, we model a range of sill velocities (4.5-6.0 km s −1 ) to estimate possible intrusion thicknesses (Smallwood and Maresh, 2002). Velocity data in the water column and the shallowest sedimentary strata were not recorded, so we assume values of 1.5 km s −1 between 0 and 64 m (i.e., seawater) and 1.8 km s −1 between 64 and 385 m (i.e., near-seabed sediments) (Figure 4).

Petrological Description of the Resolution Sill
The petrological description of the Resolution Sill was provided by Dr. G. A. Challis of the New Zealand Geological Survey (Milne et al., 1975). Based on 5 m-spaced cuttings collected between 1,911.5 and 1,958 m, the Resolution Sill is best described as a medium-to-coarse grained quartz gabbro comprising plagioclase, quartz, titanaugite, and aegerine. Minor amounts of magnetite, ilmenite, and biotite also occur. Some fine-grained, glassy, black rock chips, which contain white spherules, originate from the top contact chilled margin (see below).
At the top of the continuous core collected from the Resolution Sill, which corresponds to a depth of 1,958.2 m, the intrusion is a coarse-grained quartz syenogabbro primarily comprising titanaugite rimmed by aegerine augite, zoned plagioclase (labradorite to oligoclase) rimmed by anorthoclase, and ilmenite; fine-grained, quartz, biotite, apatite, and chlorite also occur ( Table 1). Below 1,958.3 m, quartz is absent and the Resolution Sill can be broadly classified as a teschenite that consists of plagioclase, titanaugite, analcite, anorthoclase, and occasionally olivine with accessory apatite, ilmenite, magnetite, and zeolites (Table 1). Variations in the abundance of olivine and titanaugite between ∼1,958 and 1,963 m indicate the Resolution Sill is subtly layered ( Table 1).
Frontiers in Earth Science | www.frontiersin.org

Seismic Reflection Data
This study utilizes three, zero-phase, time-migrated, 2D seismic reflection surveys (the ANZ, CB82, and Sight surveys; Figure 2). We focus on an area that covers ∼3,000 km 2 and has a total seismic line length of ∼484 km (Figure 2). Line spacing for the different vintage seismic data ranges from 3.5 to 16 km ( Figure 2B). Seismic data are displayed with a zero-phase SEG normal polarity; a downward increase in acoustic impedance correlates to a positive (black/red) reflection, whilst a negative (white/blue) reflection corresponds to a downward decrease in acoustic impedance ( Figure 5). Interval velocities derived from borehole data were used to convert the seismic reflection data from depth in seconds two-way time (TWT) to depth in meters (Figures 4, 5). We only depth-converted data above Top Basement because the lithology and physical properties (e.g., V p ) of the underlying Torlesse Supergroup are unknown ( Figure 5C).

Data Resolution
The resolution of a studied interval in seismic reflection data is dependent on the dominant wavelength (λ) of the seismic Frontiers in Earth Science | www.frontiersin.org waves, with λ = v/f, where v is he interval velocity and f is the dominant frequency (Brown, 2011). In order to distinguish reflections emanating from two distinct boundaries (e.g., the top and base of a sill), their vertical distance needs to exceed the limit of separability (∼λ/4) for the data (Brown, 2011). If the vertical distance between the boundaries is less than the limit of separability, the two reflections will interfere on their return to the surface and cannot be deconvolved; they will appear as tuned reflection packages, the true thickness of which cannot be determined (Brown, 2011). The limit of visibility (∼λ/30) defines the minimum vertical distance between two boundaries required to produce a tuned reflection package that can be distinguished from noise in the seismic reflection data (Brown, 2011). Interval velocities of 2.2-3.2 km s −1 for the sedimentary sequence in the section of interest (Figure 4), coupled with a seismic dominant frequency that decreases with depth from ∼40 to 25 Hz, suggests that the limits of separability and visibility for the data decrease with depth from ∼32 to 14 m and ∼4 to 2 m, respectively. Assuming the entire Resolution Sill has an interval of velocity of 5.2 km s −1 , equal to that of the upper 50 m intersected by Resolution-1, a surrounding dominant frequency of ∼25 Hz indicates its limits of separability and visibility are ∼52 and ∼7 m, respectively. However, if we consider that the average interval velocity of the Resolution Sill is more variable (i.e., 4.5-6.0 km s −1 ), the maximum limits of separability and visibility may be ∼60 and ∼6 m, respectively. Reflections from the top and base of the Resolution Sill where it is >60 m thick will therefore be distinguishable in the seismic data, whereas parts of the sill <60 m thick but over >6 m thick will be expressed as tuned reflection packages (see Smallwood and Maresh, 2002;Magee et al., 2015;Eide et al., 2018). Where the Resolution Sill is <7 m thick, it is unlikely to be detectable in the seismic reflection data.

Seismic Interpretation
The study area contains several high-amplitude reflections that are laterally discontinuous and can typically be subdivided into a strata-concordant inner region surrounded by a transgressive, inward-dipping limb; i.e., they display a saucershaped morphology (Figure 5). We mapped these reflections and interpret them as sills because: (i) one corresponds to the Resolution Sill intersected by the Resolution-1 borehole ( Figure 5); and (ii) they are geometrically similar to igneous saucer-shaped sills observed in the field and imaged in other seismic reflection datasets (e.g., Thomson and Hutton, 2004;Planke et al., 2005;Polteau et al., 2008a;Magee et al., 2016a). For all sills we mapped the top contact (TS) and, where seismically resolved, the base sill (BS) (Figure 5). In addition to sills, we mapped nine key seismic horizons and tied them to the Resolution-1 borehole (Figures 3-5 The limited resolution of the seismic reflection data means we cannot ascertain whether erosion has modified the geometry of the fold top (i.e., H6) and reduced its amplitude postemplacement (e.g., Figure 5C). We therefore measure amplitude along the prominent intra-fold horizon H3 (e.g., Figure 5C).
To determine fold amplitude we measure the vertical distance between the top of H3 and an inferred pre-fold datum constructed by extrapolating the regional trend of H3 from areas where there are no sills or forced folds (Figure 5C inset). The maximum vertical distance between H3 and the pre-fold datum is the maximum fold amplitude (F max ; Figure 5C). Sill thickness is measured as the vertical distance between TS and BS, with the maximum sill thickness defined by T max (Figure 5C).

Decompaction and Backstripping
Whilst several studies suggest cases where F max /T max << 1 reflects magma accommodated by elastic and inelastic deformation processes, they do not quantitatively evaluate the role of burial and compaction in modifying forced fold geometry (Jackson et al., 2013;Magee et al., 2013aMagee et al., , 2017. Loading of sedimentary sequences during burial promotes progressive loss of porosity with depth (i.e., compaction), and causes beds to become thinner and structures (e.g., faults) to flatten. The compaction of strata at any given depth is controlled by its lithology and lithostatic load. Because crystalline intrusive rocks have virtually no porosity and can be considered incompressible, T max will not change with burial. However, compaction of the overlying sedimentary sequence is expected to reduce F max and therefore decrease of F max /T max . The sedimentary sequence adjacent to the sill is overlain by a thicker column of sediment/rock, meaning it will compact more than where it is folded above the sill; this variation in lithostatic load across the fold can promote differential compaction (Hansen and Cartwright, 2006;Schmiedel et al., 2017). Evaluating the role of post-emplacement compaction in modifying forced folds is critical to establishing the relationship between the original maximum fold amplitude (F0 max ) and intrusion thickness, which can be used to inform interpretation of emplacement mechanics. To extract F0 max , we decompact and backstrip the forced fold. Note we do not take into account processes that may alter sill thickness (e.g., contraction during crystallization; Caricchi et al., 2014) and thus assume T max = T0 max .

Forward Modeling
Decompacting and backstripping sedimentary sequences imaged in depth-converted seismic reflection data involves restoring the initial porosity (ϕ 0 ) of strata at the top of the sequence from its current porosity (ϕ), by removing its overburden. This technique normally involves estimating a porosity log from sonic log data using either the Wyllie time-average method or Raymer-Hunt-Gardner empirical relationship (Wyllie et al., 1956;Raymer et al., 1980). However, given the shallow depth of our interval of interest (i.e., 1-2 km) and the limited log data available (e.g., there is no density log), we cannot reliably assess the accuracy of current porosity logs derived from these methods. We therefore apply forward modeling techniques to establish whether plausible decompacted and backstripped scenarios are realistic and how they impact fold geometry. In particular, based on the lithological information from Resolution-1, we model a series of different parameter ranges and combinations to assess potential variations between sill thickness and the original fold amplitude. Because estimates of ϕ 0 and the compaction length scale (λ), which is the inverse of the compaction coefficient, are not available, we model a range of realistic values (Sclater and Christie, 1980): (i) ϕ 0 is considered to range from 0.7 to 0.25, consistent with a range of siliciclastic sequences; and (ii) λ ranges from 3.7 to 1.4 km.
Forward modeling involved the standard back-stripping procedure described by Sclater and Christie (1980). The most common function used to model porosity decay with depth is the exponential relationship: where z is depth below seafloor. Considering the rock matrix fraction (m) complements porosity (Smallwood, 2009), whereby we can relate the pre-and post-compactional stratal thicknesses via conservation of mass to give: where z1 and z2 are the initial top and base depths of the sediment package in question (i.e., the present day depth, below the seafloor, of the top forced fold and top sill, respectively), and z3 and z4 are the original depth below the seafloor of the top and base of the sedimentary package (i.e., the original depth of the forced fold is at the seafloor, hence z3 = 0, and z4 is the depth of the sill we solve to find). Substituting equations 1 and 2 into equation 3 and integrating gives: We solve equation 4 numerically to find the thickness of the original thickness of the folded sedimentary section (i.e., z4 -z3) along individual vertical traces; by calculating equation 4 along vertical traces across the width of the sill-fold pair, we recover the original geometry of the forced fold. Because we do not know the correct input values for ϕ 0 or λ, we calculated multiple iterations of equation 4 using different, realistic combinations of ϕ 0 (0.7-0.25) and λ (3.7-1.4). This method assumes the folded layer has no elastic thickness, which we consider reasonable because much of the folded section was likely unlithified during deformation. Assuming forced folds have very little, or no, elastic thickness is consistent with forced folds having similar diameters to underlying intrusions (e.g., Hansen and Cartwright, 2006;Jackson et al., 2013).

Resolution Sill Well-Log Response
The Resolution Sill is characterized by an abrupt increase in V p , from ∼3.0 km s −1 in the overlying strata to ∼5.2 km s −1 within the sill (Figure 4). Within the sill itself, values of V p , GR, and SP vary substantially on a meter to decameter-scale (Figure 4).

Geometry
The Resolution Sill is observed on two seismic lines, with its top corresponding to a high-amplitude, positive polarity reflection (TS; Figures 4, 5). Where the base of the sill is resolved, it is characterized by a discrete, moderate-to-high amplitude, negative polarity reflection (BS) that appears to coincide with the top of the basement (TB) at a present day depth of ∼2 km (Figure 5). Although the Resolution Sill can only be mapped on two 2D seismic lines, the constraints these data provide on lateral sill tip locations allow us to interpolate its 3D geometry within a sill outline derived from comparison to exposed and seismically imaged sills (e.g., Chevallier and Woodford, 1999;Planke et al., 2005;Hansen et al., 2008;Polteau et al., 2008b;Magee et al., 2014Magee et al., , 2016a. Overall, the 54 km 2 sill has an elliptical, saucer-shaped morphology with a NW-trending, long axis of ∼6.2 km and a NEtrending short axis of ∼2.8 km (Figure 6). The strata-concordant inner sill is sub-circular, with a diameter (D) of ∼2.2 km, passing laterally into gently (8 • ), inward-dipping, up to ∼0.4 km high transgressive limbs to the SE and NW (Figures 5, 6). Toward the south-eastern edge of the Resolution Sill, at its shallowest level, the transgressive limb transitions into a strata-concordant outer rim (Figures 5, 6). Intrusion thickness appears variable across the strataconcordant inner sill, although there is a first-order decrease away from the center; assuming an interval velocity of 5.2 km s −1 for the entire sill, its thickness ranges from ∼138 m (T max ) to ≤52 m (Figures 6C, 7). Superimposed onto this outwardthinning trend within the inner sill are apparently several abrupt changes in sill thickness (e.g., there is a ∼75 m change at A-A ′ ; Figure 7). However, because the lower portion of the sill is not intersected by Resolution-1, we do not know if it is characterized by similar velocities. We also do not know whether the interval velocity of the sill varies laterally. We therefore calculate sill thickness using a range of feasible interval velocities (i.e., 4.5-6.0 km s −1 ). The envelope calculated for this velocity range constrains how thickness may vary along-strike when the sill velocity across (i.e., vertically and laterally) the intrusion is constant (e.g., 4.5 or 5.2 km s −1 ) or variable (e.g., if the velocity decreases toward its edges) (Figure 7). For a range of interval velocities, we show the sill: (i) could be up to ∼187 m thick (i.e., T max ); (ii) maintains a first-order decrease in thickness from its center outwards; and (iii) thickness still appears to show local, abrupt variations, although the magnitude of these changes may be limited depending on how velocity varies laterally (Figure 7). For example, dependent on the sills velocity configuration, the thickness change at A-A ′ could be up to ∼149 m, or down to ∼17 m. The outer portions of the transgressive sill limbs are defined by tuned reflection packages, such that their vertical thickness cannot be measured; where tuning occurs we consider intrusion thickness can range from 60 to 6 m (i.e., the limits of separability and visibility, respectively) (Figures 5, 7).
The Resolution Sill is bordered to the SW, NW, and NE by three additional saucer-shaped sills; we are unable to constrain the 3D geometry of these neighboring sills because they cannot be  Figure 5. The selected sill outline is constrained by the seismic reflection data and assumes the sill likely has an elliptical shape, similar to sills observed elsewhere (see Hansen et al., 2008). (C) Thickness map of TS-BS, i.e., where both horizons can be seismically resolved, assuming a constant sill interval velocity of 5.02 km s −1 .
mapped sufficiently on multiple seismic lines (Figure 5). In crosssection, the sills to the SW and NW display similar geometries and emplacement depths to the Resolution Sill, whereas the base of the north-eastern sill broadly coincides with Horizon H1 (Figure 5).

Host Rock Structure
Strata directly above the Resolution Sill, up to H6, are folded (Figures 5, 8). The ∼58 km 2 , elliptical (i.e., ∼6.2 km × 3 km) dome is relatively flat-topped, with uplift primarily accommodated by monoclinal bending directly above the transgressive limbs of the Resolution Sill, which cross-cut the lowermost folded strata (Figures 5, 8). The top of the fold corresponds to H6, i.e., the ∼12.5 Myr old base of the 9 m thick Harper Hills Basalt, and is onlapped by overlying, sub-horizontal strata of the Tokama Siltstone (Figures 3, 5). Whilst these seismic-stratigraphic onlap relationships indicate H6 represented the syn-intrusion free surface, the limited resolution of the seismic reflection data means we cannot ascertain whether erosion has subtly modified the geometry of the fold crest. The maximum fold amplitude (F max ) at H3 is ∼127 m, with amplitude gradually and smoothly decreasing toward the fold periphery (Figure 7). The vertical distance between H6 and TS is ∼0.75 km ( Figure 5C). Similar folds are developed above the three sills neighboring the Resolution Sill; the top of these folds all occur at H6 and their boundaries directly overlie lateral sill tips (Figures 5, 8). Although these forced folds merge in places to form compound forced folds (see Magee et al., 2014), depth-conversion of available 2D data suggests our inferred pre-fold datum captures the regional trend of the folded strata (e.g., Figure 5C). The supra-sill fold to the SW of the Resolution Sill is associated with several mound-like structures marked by moderate-amplitude, positive polarity (black) reflections that downlap onto Horizon H6 and themselves are onlapped H6-H7 strata ( Figure 5A). These mounds are up to ∼315 ms TWT high (their height in meters cannot be calculated without knowledge of their V p ) and have diameters up to ∼3.5 km. The mounds appear to have erosional bases that truncate underlying strata, including H6 ( Figure 5A). Internal reflections within the mounds are relatively poorly imaged but appear to have a convex-upwards morphology ( Figure 5A).

Fold Amplitude Compared to Sill Thickness
The maximum sill thickness (T max ) is estimated to be ∼138 m, but may range from ∼117-187 m thick depending on the interval velocity of the entire sill (Figure 7). The maximum fold amplitude (F max ) measured at H3 is ∼127 m (Figure 7). T max may thus be up to 47% greater than F max . Comparing these intrusion and fold measurements suggests F max /T max is ∼0.92, potentially ranging from ∼0.68-1.09. We also note there is a lateral offset of ∼400 m between the locations of F max and T max (Figure 7). Fold amplitude and sill thickness both display a first-order decrease toward their peripheries, although sill thickness does appear to vary abruptly in places where fold amplitude does not (Figure 7). It is difficult to determine how fold amplitude relates to the thickness of the transgressive limbs because the latter are only expressed as tuned reflection packages so only their maximum (i.e., the limit of separability, 60 m) and minimum (i.e., the limit of visibility, 6 m) thicknesses can be constrained (Figures 5, 7).

Decompaction and Backstripping Results
Decompaction of the amplitude profile across the top of the fold intersected by Resolution-1 reveals that its shape is maintained but its maximum amplitude increases from ∼127 m (i.e., F max ) to up to ∼131-185 m (i.e., F0 max ; Figure 9). Uncertainties in the decompaction input parameters means the original fold amplitude profile cannot be absolutely determined. Although calculated F0 max /T max values range from 0.70 to 1.58, the breadth of which is a function of the broader range of possible scenarios compared to F max /T max , it is clear there is a greater overlap between likely sill thicknesses and amplitude values after decompaction (Figure 9). Following decompaction, the vertical distance between H6 and TS (i.e., the emplacement depth) is ∼0.8 km.

Timing of Sill Emplacement and Forced Folding
The top of the forced fold overlying the Resolution Sill corresponds to Horizon H6, where a thin tuff, which is genetically related to the Harper Hills Basalt, is interbedded with the Tokama Siltstone ( Figure 5). Onlap of the marine, middle-to-late Miocene Tokama Siltstone onto Horizon H6 suggests it formed the paleo-seabed during sill emplacement and forced folding. Where exposed onshore, the tholeiitic Harper Hills Basalts have K-Ar ages ranging from 13.5 ± 0.4 to 11.5 ± 0.3 Ma (Sewell and Gibson, 1988), which can be used as a proxy for the age of H6. This potential age range for H6 overlaps with the radiometric date obtained for the Resolution Sill (i.e., 12 ± 2 Ma; Milne et al., 1975), suggesting sill emplacement and forced folding occurred ∼12 Ma.
Sills and forced folds adjacent to the Resolution Sill display similar seismic-stratigraphic relationships (i.e., H6 marks the fold tops) and, in places, are overlain by mound-like features we interpret as volcanoes based on: (i) their moderate-tohigh amplitude, positive polarity top contacts indicative of a downward increase in seismic velocity and density, consistent with a transition from sedimentary to igneous rocks (e.g., Symonds et al., 1998;Planke et al., 2005); (ii) observed truncation underlying strata, similar to eye-shaped hydrothermal vents, suggesting they formed via explosive activity (e.g., Jamtveit et al., 2004;Planke et al., 2005;Hansen, 2006;Magee et al., 2016b); and (iii) they have similar geometries and internal architectures to volcanic vents and volcanoes observed in other sedimentary basins (e.g., Symonds et al., 1998;Jackson, 2012;Magee et al., 2013b). Overall, our seismic-stratigraphic observations, coupled with radiometric dating of the Resolution Sill and the Harper Hills Basalt onshore, indicate a phase of magmatism and volcanic activity across the northern Canterbury Basin during the Mid-Miocene (Sewell and Gibson, 1988;Timm et al., 2010).

Emplacement Mechanics and Burial-Related Compaction
For shallow-level sills and laccoliths accommodated purely by elastic bending of the overburden, we may expect the original fold amplitude, measured at the fold top, and sill thickness to be broadly equal (i.e., F0 max /T0 max = 1) (e.g., Pollard and Johnson, 1973;Fialko and Simons, 2001;Hansen and Cartwright, 2006;Jackson et al., 2013). The ratio between the original fold amplitude and sill thickness is partially controlled by the ratio of the inner sill diameter (D) and depth of emplacement (d), with larger sills intruded at shallower depths capable of generating more bending, and thus uplift of the contemporaneous free surface, than a smaller sill at greater depths (e.g., Pollard and Johnson, 1973;Fialko and Simons, 2001;Hansen and Cartwright, 2006;Jackson et al., 2013). In particular, if the D/d ratio is >4, it is considered that the overburden will not resist bending and elastic deformation will therefore fully accommodate magma emplacement Hansen and Cartwright, 2006). Decompaction of our data indicates the Resolution Sill, which has an inner sill diameter (D) of ∼2.2 km, was emplaced at a depth (d) of ∼0.8 km beneath the contemporaneous surface (i.e., H6) and thus had a D/d ratio of ∼2.75. Given a D/d ratio < 4, the model of Pollard and Johnson (1973) suggests that the overburden may have resisted bending and, in addition to elastic deformation, promoted inelastic vertical compaction or ductile strain, thereby suppressing forced fold amplitude (i.e., F0 max < T0 max ) (see also Hansen and Cartwright, 2006;Jackson et al., 2013).
In contrast to previous studies, we quantitatively assess the impact post-emplacement compaction during burial has on fold geometry (principally amplitude) and, therefore, F0 max /T max as opposed to F max /T max (cf. Hansen and Cartwright, 2006;Jackson et al., 2013;Magee et al., 2013aMagee et al., , 2017. We show that following decompaction and backstripping, the overall fold geometry is maintained but its amplitude increases from 127 m (F max ) up to 131-185 m (F0 max ) (Figure 9). These potential F0 max values, coupled with a sill thickness of 117-187 m, means F0 max /T max ranges from 0.70 to 1.58; this is greater than our measured F max /T max range (i.e., 0.68-1.09), which we attribute to the broader range of scenarios we test in calculating F0 max /T max . Considering uncertainties in the various parameters controlling F0 max /T max measurement (e.g., sill and strata interval velocities, incorrect extrapolation of the pre-fold datum), our calculated F0 max /T max range of 0.70-1.58 suggests: (i) fold amplitude and sill thickness could be equal; (ii) fold amplitude may be less than sill thickness by up to ∼30%, a scenario consistent with a D/d ratio of ∼2.75; or (iii) fold amplitude is greater than sill thickness by up to ∼37%, which could occur when thin (i.e., with FIGURE 7 | Plot of amplitude across the fold, at H3, measured directly from the seismic reflection data (i.e., Figure 5C). We also show a range of sill thicknesses, for different seismic interval velocities, across the intrusion where TS and BS can be distinguished; a sill thickness profile considering a seismic interval velocity of 5.2 km s −1 is particularly highlighted because this is the average interval velocity for the upper 50 m of the intrusion where it is intersected by the Resolution-1 borehole. Thicknesses are not shown where the sill corresponds to a tuned reflection package along the inclined limbs, but we do highlight the maximum (max.) limit of separability and minimum (min.) limit of visibility for the sill. Note the lateral offset of F max and T max . thicknesses below the limit of visibility) sills that contributed to uplift are not seismically resolved (Reeves et al., 2018). Although uncertainties mean we cannot ascertain the true, original sillfold relationships, there qualitatively appears to be a better fit between the potential ranges of sill thickness and decompacted fold amplitude (Figure 9).
In addition to burial-related compaction, it is also worth highlighting that F max /T max discrepancies could be attributed to (Hansen and Cartwright, 2006;Jackson et al., 2013;Magee et al., 2013a): (i) reduction of fold amplitude due to erosion of the fold crest; (ii) incorrect depth conversion; (iii) strain interference with adjacent folds during deformation; (iv) out-ofplane deformation; or (v) changes in intrusion geometry. We measure amplitude from an intra-fold horizon (i.e., H3), so discount erosion of the fold crest as a mechanism for reducing F max and producing F max /T max ratios <1 (e.g., 0.81) (Figure 5C).  Figure 5C, highlighting how the measured fold shape and amplitude changes if the seismic data is decompacted and backstripped.
By using velocity data from Resolution-1 and calculating sill thickness for a range of velocity values means we have better control on depth conversion parameters than previous studies, yet our results highlight F max /T max discrepancies <1 are still plausible (cf. Hansen and Cartwright, 2006;Jackson et al., 2013;Magee et al., 2013aMagee et al., , 2017. The Resolution Sill and overlying forced fold are adjacent to and abut sill-fold pairs, of the same age, to the NW and SW (Figures 5, 8), implying strain interference between the folds may have enhanced or inhibited fold growth; however, we cannot quantify whether strain interference had a positive of negative impact on F max , particularly without access to 3D seismic reflection data. Abrupt, localized variations in thickness across the sill are not reflected by the overlying fold shape (Figures 5, 7, 9); this local decoupling between sill and fold shape may suggest vertical displacement induced by sill intrusion is distributed across an area because the overburden has some flexural strength (see Stearns, 1978). The fold profile we measure thus does not capture and may have been modified by localized out-of-plane deformation or changes in intrusion geometry; these observations imply that relatively simple ground deformation patterns may be generated by intrusions with complex geometries. Whether folded strata respond (i.e., deform) to small-scale changes in intrusion thickness is a function of emplacement depth and various host rock properties (e.g., flexural rigidity, bed thickness, co-efficient of friction between layers) (e.g., Stearns, 1978).
Overall, our work implies that explicitly accounting for burial-related compaction likely reduces measured F max /T max discrepancies (Magee et al., 2013a;cf. Jackson et al., 2013). We show that emplacement of the Resolution Sill was principally accommodated by elastic bending of the overburden but, without sufficient borehole data to accurately constrain the original fold geometry, cannot confirm whether inelastic deformation also generated space for the intrusion. Further work is required to test the impact of burial-related compaction on the geometry and amplitude of seismically imaged forced folds, and to determine how burial-related compaction signals can be deconvolved from intrusion-related inelastic deformation processes that modify fold geometry.

CONCLUSIONS
Elastic bending and uplift of overlying rock and sediment, and potentially the free surface, can accommodate emplacement of shallow-level, tabular intrusions; this intrusion-induced deformation is a form of "forced folding." Many numerical and analytical models examining sill and laccolith emplacement, as well as inversions of ground displacement data at active volcanoes used to recover information pertaining to subsurface magma movement, typically only incorporate elastic processes and neglect inelastic deformation mechanisms. Whilst the assumption that host rock deformation is purely elastic can be applied to many scenarios, several seismic reflectionbased studies have suggested that synchronous elastic and inelastic processes can generate space for magma intrusion. This interpretation that elastic and inelastic processes can accommodate magma, which is supported by some outcrop data and analytical modeling, is based on some seismically imaged forced folds having amplitudes much smaller than the thickness of the underlying intrusion; i.e., elastic bending is expected to produce folds with amplitudes broadly equal to the thickness of the underlying intrusion. However, these seismic-based studies do not quantitatively account for post-emplacement, burialrelated compaction of forced folds, which may be expected to reduce their amplitude. Through analysis of the Resolution Sill and its overlying forced fold, imaged in seismic reflection data offshore eastern New Zealand and intersected by the Resolution-1 borehole, we present the first robust decompaction and backstripping of an intrusion-induced forced fold to constrain its original geometry. Our results highlight the forced fold had an original amplitude of ∼131-185 m, but burial-related compaction has reduced its amplitude to ∼127 m. The top and base of the Resolution Sill are seismically distinguishable across its center, where it has a maximum thickness of 117-187 m, depending on the interval velocity of the entire sill. Although uncertainties still exist, we show that decompaction reduces and potentially fully accounts for apparent discrepancies between fold amplitudes and sill thicknesses. Our observations also suggest relatively simple fold shapes may be produced by complex intrusion geometries, involving local abrupt changes in thickness. Seismic reflection data provides unprecedented insights into the 3D geometry of natural intrusions and forced folds, but we highlight the need to consider the role of burial-related compaction in modifying fold shapes and amplitudes.

AUTHOR CONTRIBUTIONS
CM designed the study, carried out the seismic interpretation, interpreted the data, and wrote the manuscript. MH conducted decompaction and backstripping analysis, aided data interpretation, and edited the manuscript. CJ and SJ contributed to data interpretation and edited the manuscript.

FUNDING
Open access publication costs to be funded by University of Leeds library.