Onshore drilling by Ocean Drilling Program (ODP) Legs 150X and 174AX and offshore drilling by Integrated Ocean Drilling Program (IODP) Expedition 313 provides continuous cores and logs of seismically imaged Lower to Middle Miocene sequences. We input ages and paleodepths of these sequences into one-dimensional backstripping equations, progressively accounting for the effects of compaction, Airy loading, and thermal subsidence. The resulting difference between observed subsidence and theoretical thermal subsidence provide relative sea-level curves that reflect both global average sea level and non-thermal subsidence. In contrast with expectations, backstripping suggests that the relative sea-level maxima in proximal onshore sites were lower than correlative maxima on the shelf. This requires that the onshore New Jersey coastal plain has subsided relative to the shelf, which is consistent with models of relative epeirogeny due to subduction of the Farallon plate. These models predict subsidence of the coastal plain relative to the shelf. Although onshore and offshore sea-level estimates are offset by epeirogeny, the amplitude of million-year–scale Early to Middle Miocene sea-level changes seen at the New Jersey margin is generally 5–20 m and occasionally as great as 50 m. These events are interpreted to represent eustatic variations, because they occur on a shorter time frame than epeirogenic influences. Correction for epeirogenic effects largely reconciles differences between onshore and offshore relative sea-level estimates and suggests that backstripping provides a testable eustatic model for the Early to Middle Miocene.

One of the outstanding challenges in studying Earth history is documenting the timing and magnitude of global sea-level change (e.g., Haq et al., 1987, 1988; Miller et al., 2005). Here we use eustasy to mean the global change in sea level relative to a fixed point, e.g., the center of the Earth (Posamentier et al., 1988). We use the term relative sea-level change to include changes in eustasy coupled with epeirogenic (broad regional uplift; Grabau, 1936) and local changes in the height of lithosphere relative to the center of the Earth. That is, relative sea-level change is measured relative to a fixed point on the crust (Posamentier et al., 1988). Eustasy is of particular importance because it serves as the datum against which Earth’s tectonic and climate history can be measured (e.g., Miller et al., 2005). In this paper, we attempt to untangle thermal subsidence and epeirogeny from eustasy.

The timing and magnitude of relative sea level has been studied in detail at the Late Cretaceous to Miocene of the New Jersey margin. Onshore New Jersey drilling by Ocean Drilling Program (ODP) Legs 150X and 174AX in concert with Legs 150 and 174A outer shelf and continental slope drilling (Fig. 1; e.g., Miller et al., 1997, 1998a; Mountain et al., 2010) has shown that the timing of sequence boundaries (erosional surfaces recognized in cores and seismic profiles), is consistent with δ18O increases from deep-sea records (Miller et al., 1991, 1996a, 2005, 2011; Browning et al., 2008). This suggests that relative sea-level changes observed on the New Jersey margin were caused, at least in part, by eustatic variations due to ice growth and decay.

Backstripping is a modeling technique that accounts for the effects of sediment compaction, sediment loading, and, in this case, thermal subsidence. Applied to the onshore coreholes, backstripping provides a relative sea-level record for the New Jersey margin that in the absence of other tectonic effects yields a testable estimate of eustatic change (Miller et al., 2005; Kominz et al., 2008). However, mantle tomographic studies coupled with models of lithospheric epeirogeny suggest that the New Jersey margin has undergone broad tectonic subsidence over the past 50 million years (e.g., Conrad et al., 2004; Moucha et al., 2008; Spasojević et al., 2008). Thus, more work is required to separate eustatic and epeirogenic effects in this region.

Several processes must be accounted for to untangle eustasy from epeirogeny. Most backstripping estimates of the magnitude of sea-level change in this region have used a one-dimensional approach that assumes an Airy response to sediment loads (Kominz et al., 1998; Miller et al., 1998a, 2005; Van Sickel et al., 2004; Kominz et al., 2008). Only one model, of Upper Eocene to lowermost Miocene strata, used a two-dimensional backstripping approach that incorporates flexural rigidity of the lithosphere to account for subsidence caused by sediment loads at a distance (Kominz and Pekar, 2001). There are several complications inherent in estimates of New Jersey margin relative sea level. One issue is the fact that coastal plain sediments rarely contain a complete record of sea-level change. They generally preserve only the transgressive and highstand systems tracts, leaving lowstand sediments farther seaward beneath what is now the continental shelf (Miller et al., 1998a). As discussed above, other tectonic effects have been postulated in this region so that tectonic subsidence may not be entirely thermal. In particular, the arrival of the Farallon slab beneath the North American east coast ca. 75 Ma requires that New Jersey subsided beyond the predicted thermal subsidence (Conrad et al., 2004; Kominz et al., 2008; Moucha et al., 2008; Müller et al., 2008; Spasojević et al., 2008). This means that the stratigraphic succession in this region has been imprinted by eustatic, thermal, and epeirogenic processes. Additionally, the whole Earth response to loading effects of water coupled with the unloading of glaciers (glacial isostatic adjustment [GIA]) has been shown to vary globally (e.g., Peltier, 1998) with impact in this region (e.g., Raymo et al., 2011). This effect is most pronounced during the large Northern Hemisphere ice ages of the past 2.7 m.y., but GIA influences the reference frame of older records as well (e.g., Raymo et al., 2011).

Drilling data provided by IODP Expedition 313 (hereafter “Exp 313”) on the New Jersey shallow continental shelf focused on Lower to Middle Miocene strata (Mountain et al., 2010). This data set presents an opportunity to estimate amplitudes of offshore relative sea-level change based on cores, logs, and seismic profiles that can be compared to the onshore results. By providing estimates of the magnitude of relative sea-level change, we can begin to address the magnitude and timing of both eustasy and more regional epeirogeny.

Exp 313 drilled three coreholes (sites M27, M28, and M29; Fig. 1) in ∼30 m of water, targeting Miocene sequences that also were cored in multiple locations on the onshore coastal plain. The Miocene section is well imaged on a grid of seismic profiles that show a series of clinothems, which are packages of sediment that prograde seaward and are bounded by surfaces (in this case, sequence boundaries) with distinct sigmoidal (clinoform) geometries (Fig. 2; Mountain et al., 2010). These sites, drilled as a transect along seismic line Oc270 529, provide a two-dimensional cross section of several stratigraphic sequences formed between 12 Ma and 22 Ma (Mountain et al., 2010). Here, we present the results of a one-dimensional backstripping study of this new data set. As such, it is directly comparable to the one-dimensional modeling already published from the coastal plain (Kominz et al., 2008). By extending that study to this offshore location, it is possible to provide a preliminary estimate of the magnitude of million-year–scale relative sea-level changes and to consider the effects of the Farallon slab and GIA on relative sea-level change across the New Jersey margin.

Backstripping utilizes compaction, age, and paleodepth observations from cores or outcrops to estimate how the basement would have subsided (tectonic subsidence [TS]) in the absence of sediments and eustatic sea-level change (ΔSL) (Steckler and Watts, 1978; Bond and Kominz, 1984):
where: ρ is density of the asthenospheric mantle (m), the decompacted sediment (s), and seawater (w); WD is local water depth; S* is the decompacted sediment thickness; and φ is the basement response function.
In practice, backstripping begins with measuring porosity as a function of depth and lithology in a borehole or outcrop. These data are used to estimate the porosity, and thus, decompacted thickness (S*) of every part of the sedimentary column through time. The sediment column is removed at each time step and replaced by a column of seawater that represents the hole that sediment of the estimated porosity, thickness, and density would have filled. In calculating the effect of the sediment load, we have applied a one-dimensional approach that assumes the underlying lithosphere has no rigidity (the Airy model, in which φ, in Equation 1 is taken as 1). We also account for changes in water depth as described below. For that, paleodepth estimates (WD) are added to the water column that was filled with sediment (above) to determine the total subsidence of the basement in the absence of sediment. The result is termed the first reduction, or R1 (Bond et al., 1989):
where the variables are the same as described in Equation 1. R1 does not take into account global average sea-level change ( = eustasy; ΔSL in Equation 1). R1 is based on observed data that can be obtained from a corehole with a reasonable degree of accuracy.
Tectonism at a passive margin may be assumed to follow the theoretical behavior of a cooling plate due to lithospheric stretching during rifting, as demonstrated by the fact that long-term (10–100 m.y.) passive margin records can be largely explained with an exponential fit (e.g., McKenzie, 1978; Royden and Keen, 1980). Because of the predictable nature of tectonics in this setting, we can estimate eustasy by separating the long-term, thermal component of subsidence from any perturbations observed in the R1 curve. A second parameter, R2, is calculated as the difference between R1 and a cooling plate model fit to the observed R1 data. We fit a thermal plate model (following Kominz et al., 2008) to the R1 curve to determine the thermal component of subsidence. The difference between the estimated thermal subsidence (taken as TS, Equation 1) and R1 (Equation 2), after being corrected for the change in water load (again assuming a one-dimensional Airy response of the underlying lithosphere) yields the second reduction R2:

Employing Equations 1 and 2 to define TS and R1 explicitly defines R2 as a eustatic estimate. However, this is only true if all tectonics at the borehole in question are governed by thermal subsidence. Additional errors may arise as a result of assuming Airy isostasy, in a basin that actually responded flexurally to the sediment load because subsidence is underestimated at the locus of loading and overestimated at the periphery of the basin. However, in our approach we are looking at the difference between R1 and the theoretical thermal subsidence fit to the observed subsidence, R1. The form of subsidence resulting from the flexural response to loading of a thermal basin is thermal in both the center of the basin and its periphery (Bond et al., 1988). Therefore, the difference between R1 and thermal tectonics, R2, is not affected by one-dimensional backstripping (Kominz et al., 1998). What is lost through one-dimensional analysis of a flexurally subsiding thermal basin is the detailed relationship between data points along a two-dimensional profile that could add insights for interpretations (e.g., Steckler et al., 1999; Kominz and Pekar, 2001). Observed relative sea level is also complicated by glacial isostatic adjustment (GIA; Peltier, 1998), by gravitational, rotational, and flexural effects due to changing ice sheets (collectively known as “static equilibrium” effects; Kopp, et al., 2010) and by dynamic topography (e.g., Milne et al., 2009). These issues will be considered in evaluating the results.

Lithology and Porosity

We use the backstripping approach of Bond et al. (1989) in which porosity is assumed to be lithology dependent and the porosity of mixed lithologies is calculated based on the proportion of lithologies present. This has been shown to be a viable approach for marine sediments (Kominz et al., 2011). Porosities based on moisture and density measurements from the three sites used in this study were combined with smear slide interpretations, grain-size analysis of discrete samples, and downhole log data (Mountain, et al., 2010; Miller et al., 2013a; Ando et al., 2014) to obtain lithology-dependent porosity versus depth plots for the New Jersey shelf (Fig. 3 and Table 1). Lithologies were extended from the depths of each discrete measurement to the nearest stratigraphic boundary (Fig. 2) above and below that measurement; otherwise, the lithologic boundaries were placed between samples. Siliciclastic sand-, silt- and clay-dominated intervals are sufficiently distinct in the Exp 313 cores to provide valid comparisons with standard, lithology-based porosity versus depth plots published elsewhere. While the Exp 313 data show no consistent pattern with depth, their trends are close to the >90% clay, the >90% silt, and the >60% sand curves based on a global Ocean Drilling Program (ODP) database (Kominz et al., 2011). Thus, we use the Kominz et al. (2011) porosity versus depth relations to decompact the sediments in these coreholes (Fig. 3). High-end and low-end porosity versus depth curves that encompass the range of observed porosities were also applied for comparison (see Supplemental Table S11 and Supplemental Fig. S12).

Ages and Environments of Deposition

The focus of this work is the Lower–Middle Miocene section, for which considerable data are available and has been synthesized to provide input for modeling (Table 2; detailed input data are provided in Supplemental Table S23). Age estimates are based on Browning et al. (2013), who used calcareous nannofossil, dinoflagellate cyst, and diatom biostratigraphy, and numerous strontium (Sr) isotopic age estimates of calcium carbonate from mollusk shells, shell fragments, and foraminiferal tests. Ages were assigned using the time scale of Gradstein et al. (2012; GTS2012). By using Sr-isotope stratigraphy, Browning et al. (2013) overcame the challenges of poor magnetostratigraphy and poor biostratigraphy posed by coarse, clastic, nearshore sediments. They were able to generate age resolution of typically ±0.5 m.y., and in many sequences, age resolution was as good as ±0.25 m.y., where age resolution refers to uncertainties in correlation to the geologic time scale (Browning et al., 2013).

Superposition and sequence stratigraphic principles provide additional relative age constraints. Our biostratigraphic and Sr-isotopic data have been combined with seismic profiles, downhole logs, sedimentological data, and sequence stratigraphic interpretations to provide age estimates for sequence boundaries and several prominent subsequence stratal surfaces (Browning et al., 2013). In this study, dates of some sequence boundaries and maximum flooding surfaces (MFS) have been shifted slightly from previous work of Mountain et al. (2010) and Browning et al. (2013) to honor new, detailed sequence stratigraphic correlations by Miller et al. (2013b) and to reconcile minor differences between onshore and offshore locations (Table 3). In particular, sequence boundaries are assumed to represent times of non-deposition that may represent more time at one location than another, but the sediments below the unconformity are everywhere older than those above. Additionally, MFS (based on core data) are assumed to be correlative and are given a specific date within the sequence. Finally, where onshore and offshore sequences are correlative, minor shifts in ages (0–0.45 m.y.) of the offshore strata were made to align correlative sequences. All assigned dates are within the error limits presented by Browning et al. (2013).

Environments of deposition at the three well sites were based on benthic foraminiferal assemblages, planktonic foraminiferal abundances, and palynological proxies (Katz et al., 2013; McCarthy et al., 2013; Miller et al., 2013b) coupled with lithofacies interpretations (Browning et al., 2013). Generally the lithofacies environments are: foreshore (<10 m), upper shoreface (10–20 m), shoreface-offshore transition (20–30 m), and offshore (below storm wave base that we adopt as ≥30 m). Benthic foraminiferal bathymetric zones are defined as inner neritic (0–30 m), middle neritic (30–100 m), and outer neritic (100–200 m) (van Morkhoven et al., 1986). Miller et al. (1997) and Katz et al. (2013) used a subdivision of these zones on the New Jersey margin when they interpreted Elphidium-dominated biofacies as <10 m, Hanzawaia hughesi-dominated biofacies as 10–25 m, Pseudononion pizarrensis–dominated biofacies as 25–50 m, Bulimina gracilis–dominated biofacies as 50–80 m, and Uvigerina-dominated biofacies as >75 m. For deeper biofacies not recovered by Miller et al. (1997), Katz et al. (2013) used key taxa (e.g., Cibicidoides pachyderma, Cibicidoides primulus, Hanzawaia mantaensis, and Oridorsalis umbonatus) often found in high-diversity, low-dominance assemblages that indicate outer neritic paleodepths (100–200 m) (e.g., Parker, 1948; Poag 1981; van Morkhoven et al., 1986; Katz et al., 2003). Because changes from the Uvigerina spp. biofacies to the deeper biofacies occur gradually within sequences, Katz et al. (2013) inferred that the maximum water depths are in the shallower part of the outer neritic zone (∼120 m). Water-depth ranges as well as best estimate water depths are used in modeling.

Stratigraphy above and below the Lower–Middle Miocene

Sediments beneath the Lower to Middle Miocene section sampled by Exp 313 have been compacted under the load of the Miocene and younger sediments. This compaction must be taken into account to properly estimate R1 for the Miocene strata. Based on a contour map of depth to basement (Benson, 1984), we estimate that ∼8 km of pre-Miocene sediments were deposited beneath the three offshore sites. For modeling purposes, we use the along-strike Anchor Gas, Dickinson No. 1 rotary well (Poag, 1985; Olsson et al., 1988; Sugarman et al., 2011; Fig. 1) to estimate the lithology and age data of the underlying strata.

The Lower to Middle Miocene strata also compacted beneath Upper Miocene and younger sediments. These sediments were either poorly sampled (Exp 313 sites M27 and M29) or not sampled at all (site M28). In our model, porosity (and thus compaction) depends only on depth of burial; thus the detailed lithology and water depth of this younger part of the section have no impact on R1 results. Therefore, the decompaction of the Lower–Middle Miocene focus of the paper is valid because the thickness of overlying strata is known and included in the model. The thermal curve is fit to the R1 water-depth values based on the pre-Miocene Anchor Dickinson sediments overlain by sites M27, M28, or M29 detailed stratigraphic data that are, in turn, overlain by 209, 243, or 342 m of younger strata observed at the latter three sites, respectively. The effects of the younger subsidence history on the magnitudes of Miocene relative sea level are minor, because the R2 results are registered to present sea level, which is unaffected by the intervening upper Miocene and younger section. Corehole water depths at Sites M27, M28, and M29 are 34, 35, and 36 m below present sea level, respectively (Mountain et al., 2010). The registry to present sea level is a reasonable first assumption; although recent GIA effects on the order of 20 m may complicate absolute sea-level estimates (Raymo et al., 2011) and are taken into account in our discussion in the Results section.

Coastal Plain Relative Sea-Level Curve

In this paper, the time scale used by Kominz et al. (2008) to derive a sea-level curve from Miocene sequences has been updated to GTS2012 (Gradstein et al., 2012). The offshore results are of particular interest in comparison to the coastal plain relative sea-level estimates. However, during the two decades these drill core samples have been studied, revisions in the time scale—from Berggren et al. (1995) (BKSA95) to Lourens et al. (2004) and from the GTS2004 to the GTS2012 (Gradstein et al., 2012)—have resulted in small but significant adjustments to the Early to Middle Miocene. Early to Middle Miocene ages are largely unchanged from BKSA95 (used for onshore sequences by previous studies; e.g., Kominz et al., 2008) and offshore by Mountain et al. (2010) from GTS2004 (used by Browning et al., 2013), except for the age of the Oligocene/Miocene boundary. The Early to Middle Miocene time scale of GTS2004 differs from GTS2012 only by placement of the Langhian/Burdigalian boundary (0.16 m.y. older in GTS2012). All coreholes used in generating the coastal plain R2 curve have fairly well defined age ranges for the Miocene sequences, with age errors generally ±0.5 m.y. In sediments younger than 15 Ma, the errors can be larger (i.e., ages are constrained by Sr isotopes that have an error of ±0.75 to ±1.0 m.y.; Browning et al., 2013, their fig. 3). In downdip onshore coreholes, the age ranges are better constrained than they are at updip locations. We used the better dated sequences to constrain the less well dated coreholes. This generally resulted in age shifts of less than 0.5 m.y. and slightly greater amplitudes in the sea-level estimates, because incidences of destructive interference due to age miscorrelations have been reduced. The Cape May Zoo site (Sugarman et al., 2007) was analyzed subsequent to the Kominz et al. (2008) compilation, and it has been included in our backstripped data set (see Results section). Lithologies at all coastal plain sites were decompacted using the coastal plain porosity versus depth relations of Van Sickel et al. (2004).

Coastal Plain Relative Sea-Level Curve

Cape May Zoo R1 and R2

The Cape May Zoo corehole (Sugarman et al., 2007) provides additional Lower to Middle Miocene data that were not included in previous backstripping studies (Kominz et al., 2008). Younger sequences at this site are poorly dated, and they are not included in our R2 analysis. The Oligocene and Eocene portion of the R1 curve (Fig. 4A) is based on sediment from the Cape May Site (Miller and Snyder, 1997), with older units based on the Anchor Dickinson well (Poag, 1985; Kominz et al., 1998). The beginning of thermal subsidence is taken as 140 Ma, the time when offshore sediment loading overcame lateral thermal uplift landward of the hinge zone (Steckler et al., 1988). For the Middle Miocene section, the subsidence rate is slightly greater than that of the thermal curve, suggesting rising relative sea level during this interval and falling relative sea level between deposition of these Miocene sediments and the present (Fig. 4B).

Combined Results, Coastal Plain Coreholes

Several coreholes drilled over the past decade in the New Jersey coastal plain penetrate Lower and Middle Miocene strata, locally called the Kirkwood Formation. These include coreholes at Ancora (Miller et al., 1999), Bass River (Miller et al., 1998b), and Island Beach (Miller et al., 1994a), where only a few of the Kirkwood sequences were sampled and dated (Kw1a and Kw1b and possibly Kw2a or b) and coreholes at Ocean View (Miller et al., 2001), Millville (Sugarman et al., 2005), Cape May (Miller et al., 1996b), and Atlantic City (Miller et al., 1994b) that sampled multiple Kirkwood sequences. Early and Middle Miocene R2 curves for these sites, in addition to the new Cape May Zoo site (Sugarman et al., 2007), show that several Kirkwood sequences are well represented (Fig. 5A). In previous studies, the timing of sequences did not always match from one corehole to the next due to uncertainties in assigning ages to sequences, particularly in updip locations (e.g., the Ancora and Bass River sites). This resulted in a blurring of the sequences so that sequence boundaries were not seen as hiatuses in the combined record (e.g., Kominz et al., 2008). While it is possible that the poor match is real, in most cases the uncertainty of age dates is ∼1 m.y. and larger in updip sites. We assume that the coastal plain sequence boundaries are regionally correlated, and we accept the ages of the best dated sequences and adjust the other coreholes accordingly. Seismic correlations onshore and offshore confirm the correlations (Monteverde et al., 2008; Iscimen, 2014). In general, the Cape May Zoo sequences are well dated and provide ages for the base of most sequences. However, for Kw1 and Kw2a, the sequences at Ocean View provide the best age constraints. Sequences Kw0, Kw1c, and Kw3a are best represented at the Cape May corehole. The tuned sequence ages allow us to stack sequences from several coreholes; the composite distinguishes individual sequences in the relative sea-level record (Fig. 5B).

The composite relative sea-level curve (Fig. 5C; data are provided in Supplemental Table S34, section A) was generated following the procedures of Kominz et al. (2008). The minimum possible range of sea-level estimates is obtained as the minimum high-end sea-level estimate and the maximum low-end sea-level estimate. The best estimate relative sea level is taken as the average of the high- and low-end estimates. Note that despite age model tuning, we have not been able to separate sequences Kw1a and Kw1b that show no discernable hiatus at Ocean View or other coreholes. Nonetheless, we have separated sequences Kw3c into two sequences and separated Kw3a from Kw3b and Kw1b from Kw1c. Error ranges have been reduced in some sequences as a result of matching high and low relative sea-level curves, and the timing of the sequences has shifted slightly, mainly due to changes in the geologic time scale.

Exp 313 Offshore Sites M27, M28, and M29

A critical test of the efficacy of one-dimensional backstripping for estimating relative sea level is the degree to which the calculated magnitude of R2 is consistent among boreholes. In most instances and within errors of water-depth estimates, our study passes this test for the three Exp 313 drill sites (Fig. 6). Although using low- and high-end porosity versus depth curves increases ranges of uncertainty, the form of m.y.-scale relative sea-level variations remains unchanged (see Supplemental Fig. S25). Sequences younger than sequence boundary m4 (ca. 12.6 Ma) are included in the modeling, but they are not the focus of this study (see Table 2). Poor age and paleoenvironment control for these younger sequences could explain the lack of agreement between sites (Fig. 6). To compare these new offshore results with our relative sea-level curve generated using coastal plain boreholes, we construct a relative sea-level curve based only on the offshore sites, following the method of Kominz et al. (2008) described above. That is, the error ranges are obtained by taking the minimum high estimate and the maximum low estimate of R2. Where multiple coreholes sample the same sequence, this results in a reduction in the range of error bars, because only the overlapping magnitudes are taken to be valid (Fig. 7A; results are provided in Supplemental Table S3 (section B, [footnote 4]). The “best estimate” R2 curve is based on an average of best estimates for all coreholes that sample sediment of that age. Where R2 results from multiple coreholes do not overlap, the best estimate average does not always fall between the high and low ranges (e.g., m5.3 and m5.2; Fig. 7). In these cases, the error range has been extended to include the best estimate averages (Fig. 7B). Despite these inconsistencies, there is an overall good age agreement, with all three offshore coreholes showing similar R2 trends and magnitudes (Figs. 6 and 7). This suggests that we are reconstructing relative sea-level change at the New Jersey shallow shelf.

The relative sea-level curve generated from offshore data shows higher magnitudes of sea-level fluctuations than the onshore data (Fig. 8A). Taking into account the full error range, the largest long-term sea-level range allowable in the shelf data is 75 m (16.7–13 Ma) compared to a 55 m range in the onshore data (20.4–20.3 Ma). On the other hand, the minimum allowable sea-level ranges of the two data sets are 32 m in the offshore data (16.7–12.2 Ma) and a 17 m in the onshore data (19.2–13.4 Ma). These ranges are entirely compatible with those observed in comparable time intervals at the Marion Plateau by John et al. (2011). They report four sequences with maximum sea-level changes of 65 m and minimum ranges of 26 m between 13.8 and 16.6 Ma. All estimates are lower than the 100–120 m variations implied by the global compilation of Haq et al. (1987), but are compatible with magnitudes of 20–50 m suggested by Haq and Al-Qahtani (2005) for the Arabian Peninsula.

Comparing the averaged R2 curves from Exp 313 with the revised composite onshore relative sea-level curve has implications for models of tectonics and GIA effects (Fig. 8A). From a sequence stratigraphic approach, we would expect coastal plain sequence deposition to represent the most proximal portion of transgressive and/or highstand systems tracts. These might occur during slowly rising relative sea level that may be associated with non-deposition offshore. Thus, coastal plain sequences would be present during the maximum amplitude of relative sea level. Our initial comparison of onshore and offshore sites (onshore is shown as light gray units in Fig. 8A) shows that deposition on the coastal plain sometimes occurs during hiatuses in the corresponding offshore sedimentation sampled by Exp 313. Where this is the case, it implies that the sequence boundaries of offshore units indicate relative sea-level rise coupled with offshore sediment starvation rather than the more commonly accepted model of sea-level fall (e.g., Loutit et al., 1988). In most cases, onshore and offshore sequences show significant temporal overlap, with the offshore generally more complete. There are similar trends in offshore and onshore data in overlapping intervals. In fact, sequences Kw3b and Kw3c are, at least in part, compatible with relative sea level obtained from analyzing offshore strata (Fig. 8A). Thus, in most cases, both onshore and offshore sequence boundaries represent sea-level falls.

The R2 sea-level estimates obtained from the onshore sites (0 ± 20 m above present sea level) are considerably lower than those indicated by the offshore sites (40 ± 30 m above present sea level), particularly in the older Miocene units. This is in direct contrast with stratigraphic principles that require sea-level rise to be higher when proximal land is flooded (e.g., Posamentier and Vail, 1988). As such, it suggests that the New Jersey margin has experienced non-thermal epeirogeny, resulting in differential subsidence of the onshore and offshore portions of the margin. The dominant epeirogenic effect predicted to have occurred in this region is the overriding of the Farallon plate. As the cold, dense lithosphere of the subducted slab passes beneath the margin, coastal New Jersey experiences subsidence (Conrad et al., 2004). Moucha et al. (2008) predicted that over the past 30 m.y., the coastal plain of New Jersey has subsided 50–100 m more than offshore New Jersey. The minimum value of 50 m of differential subsidence over 30 m.y. implies a relative coastal plain uplift of 33 m at 20 Ma, reducing to 20 m at 12 Ma. We add this differential subsidence to the onshore relative sea-level curve to make it comparable to the offshore R2 results. This correction for epeirogenic effects largely reconciles differences between onshore and offshore relative sea-level estimates.

An additional variation between onshore and offshore relative sea level is predicted by the effect of ice loading (GIA). Raymo et al. (2011) estimated the impact of ice volume and GIA effects on observations of Pliocene sea-level maxima. They found that GIA effects have diminished through time for all but the last glacial maximum and the very near-field location of the melted ice sheets. However, the impact of the last glacial maximum on the elevations of the New Jersey coastal plain and the New Jersey shelf estimated by Raymo et al. (2011) will also affect older sea-level estimates by distorting today’s datum. They found that the New Jersey coastal plain yields relative sea-level estimates that are ∼12 m above eustatic sea level, while the shelf yields relative sea-level estimates that are ∼18–24 m above eustatic sea level (depending on distance offshore and on the assumptions made for lithosphere rheology). Therefore, offshore relative sea-level estimates will always appear to be 6–12 m higher than coastal plain estimates. Thus, to compare onshore relative sea-level maxima with offshore data, 6–12 m needs to be added to all onshore data. Rowley et al. (2013) also suggest five to ten meters of relative uplift occurred between the coastal plain and shelf as a result of GIA, consistent with the Raymo et al. (2011) results. The effect of the Farallon slab in conjunction with the GIA effects shift relative sea-level variations predicted by the coastal plain curve relative to the offshore R2 curve (Fig. 8B). The resulting coastal plain adjusted-R2 magnitudes, in comparison with the offshore data, are consistent with the assumption that the deposition on the onshore coastal plain generally occurs at or near the highest magnitude of sea level. Our corrected relative sea-level curves for onshore and offshore New Jersey provide a working model for eustatic changes for the Early to Middle Miocene.

It is premature to draw broad conclusions about the implications of these results for the relationship between sequences and systems tracts for sequence paradigms because of limitations in the backstripping technique and because of uncertainties in correlating sequences from onshore to offshore. Such detailed comparisons are hampered by the fact that onshore data are compiled from seven widely dispersed locations, while the offshore data set consists of three coreholes along a single dip transect.

Correlation of the data sets relies on both chronostratigraphy and the basic assumptions of sequence stratigraphy (including superposition). All strata above a sequence boundary are younger than the strata below; this is also true for maximum flooding surfaces. This allows correlation of sequences amongst onshore locations and offshore locations with finer resolution than provided by the chronostratigraphy (e.g., we are reasonably certain that sequence Kw1a is correlatable throughout the coastal plain coreholes). We are also reasonably certain of the onshore to offshore correlation not only through chronostratigraphy but also through pattern matching and seismic correlations. For example, the onshore composite sequence Kw2a sequence can be precisely correlated with composite sequence m5.4 using well logs and seismic profiles (Iscimen, 2014). Thus, though the uncertainty on the numerical age of any given backstripped estimate is large (±0.25 to ±0.5 m.y.), physical correlations (principles of superposition, seismic control, and sequence stratigraphy) mean that our relative age correlations amongst sites is finer.

If deposition did not follow the tenets of sequence stratigraphy, then correlation would rely solely on chronostratigraphy. The age uncertainty on the ages of sequences is large (±0.25–0.5 m.y. offshore and ±0.5–1.0 m.y. onshore) compared to the mean sea-level cycle duration of 0.84 ± 0.52 m.y. Exp 313 also sampled several higher-order sequences (100/405 k.y.; Browning et al., 2013); thus the average duration of our sequences is shorter than 1 m.y. Nevertheless, our sequences are of a higher order than predicted by changes in mantle dynamic topography that appear to operate on longer-term time scales (2–100 m.y.; Petersen et al., 2010). Reliance on chronostratigraphy alone would not alter the pervasively lower relative magnitude of sea level recorded on the coastal plain as compared to that observed on the shelf. That is, relative subsidence of the coastal plain would still be required, supporting the epeirogenic model. We expect more precise relative sea-level interpretations when the onshore and offshore core data are combined with seismic imaging in a full two-dimensional backstripping approach (e.g., Kominz and Pekar, 2001). Future work will include additional correlations of offshore sequences using seismic profiles and well logs (e.g., Iscimen, 2014) and two-dimensional backstripping (e.g., Steckler et al., 1999).

In summary, comparison of R2 results from the new offshore coreholes with the revised coastal plain sequences requires relative subsidence of the coastal plain since the Middle Miocene. This is compatible with modeling of mantle-driven subsidence due to the subducted Farallon slab (e.g., Moucha et al., 2008).

Backstripped relative sea-level estimates of Middle to Late Miocene sequences from the New Jersey coastal plain and three new IODP shelf sites are generally internally consistent with respect to the timing of sequence boundaries and relative sea-level variations. Coastal plain sedimentation tends to predict relative sea-level changes that are consistent with those seen offshore. Where onshore and offshore sedimentation corresponds, the rising and falling portions of the relative sea-level curve can be correlated, although the magnitude of offshore and onshore estimates is offset. This result requires that the New Jersey passive margin has undergone epeirogeny and, in particular, that the offshore shelf has subsided less than the coastal plain. Thus, it is consistent with relative epeirogeny due to subduction of the Farallon plate (Moucha et al., 2008) and due to GIA effects of the last deglaciation (Raymo et al., 2011). The amplitude of Late to Middle Miocene m.y.-scale sea-level changes seen at the New Jersey margin is generally 5–20 m but occasionally is as great as 50 m.

This manuscript was greatly improved by suggestions from one anonymous reviewer, reviewers Hugo Pouderoux and Olivier Dauteuil, and Geosphere Associate Editor Jean-Noël Proust. This work was supported by funds from the Faculty Research and Creative Activities Award, Western Michigan University to Kominz. Support was provided from National Science Foundation grants EAR-1052257, OCE-1154379, and OCE14-63759 to Miller. Funding was supplied by the U.S. Science Support Program Consortium for Ocean Leadership to Katz, Miller, Browning, and Mountain. Samples were provided by the IODP.

1 Supplemental Table S1. Porosity versus depth curves. Please visit http://dx.doi.org/10.1130/GES01241.S1 or the full-text article on www.gsapubs.org to view Supplemental Table S1.
2 Supplemental Figure S1. Integrated Ocean Drilling Program Expedition 313 porosity versus depth for sediments dominated by (>50%) sands, silts, and clay. Please visit http://dx.doi.org/10.1130/GES01241.S2 or the full-text article on www.gsapubs.org to view Supplemental Figure S1.
3 Supplemental Table S2. Ages, densities (rho), lithologies, and water depths used as input for sequences from coreholes M27 (section A), M28 (section B), and M29 (section C). Please visit http://dx.doi.org/10.1130/GES01241.S3 or the full-text article on www.gsapubs.org to view Supplemental Table S2.
4 Supplemental Table S3. Relative sea-level curve for New Jersey Shelf and relative sea-level curve for New Jersey Coastal Plain. Please visit http://dx.doi.org/10.1130/GES01241.S4 or the full-text article on www.gsapubs.org to view Supplemental Table S3.
5 Supplemental Figure S2. Relative sea-level (R2) curves from Integrated Ocean Drilling Program Expedition 313 coreholes. Please visit http://dx.doi.org/10.1130/GES01241.S5 or the full-text article on www.gsapubs.org to view Supplemental Figure S2.