2.1 Geology and glacial history of the Olympic Mountains

At present, the Juan de Fuca Plate subducts obliquely with respect to the overriding North American Plate (Fig. 1a) at 34 mm yr⁻¹ at the latitude of the Olympic Mountains (Doubrovine and Tarduno, 2008). The forearc high of the subduction zone comprises (from north to south) Vancouver Island, the Olympic Mountains, and the Oregon Coast Range, and it lies west of a forearc low (e.g., Georgia Lowlands, Puget Lowlands) and the active volcanic arc (Fig. 1a). Seismic imaging suggests a flatter subduction angle beneath southern Vancouver Island and the Olympic Mountains (Hayes et al., 2012; McCrory et al., 2012) compared to areas in the north and south (Fig. 1a). The modern configuration of the subduction zone was established by the latest Eocene (e.g., Brandon and Vance, 1992) after accretion of the Coast Range Terrane to the North American continent (Fig. 1c). This terrane represents a large oceanic plateau and extends from the southern tip of Vancouver Island to Oregon (Eddy et al., 2017; Phillips et al., 2017; Wells et al., 2014).

The accretionary wedge of the subduction zone is exposed onshore within the Olympic Mountains (Fig. 1a) and is composed of Eocene–Miocene flysch (Brandon et al., 1998; Tabor and Cady, 1978). This part of the mountain range is known as the Olympic Structural Complex (Brandon et al., 1998) and is separated from the surrounding Coast Range Terrane by the Hurricane Ridge thrust fault (HRF; Fig. 1c), a major discontinuity traceable in seismic surveys (e.g., Clowes et al., 1987; Calvert et al., 2011). Minor sedimentary rocks of Eocene age (Eddy et al., 2017; Tabor and Cady, 1978) are contained within the Coast Range Terrane besides the predominant ∼50 Ma marine and subaerial basaltic rocks (Eddy et al., 2017). Exhumation of the range commenced at 18 Ma and since 14 Ma, the orogen is presumed to be in flux steady state (Batt et al., 2001; Brandon et al., 1998; Pazzaglia and Brandon, 2001).

Plio–Pleistocene glaciation has strongly influenced the present-day appearance of the Olympic Mountains (Fig. 1b). During its maximum extent at ∼14 ka, the Cordilleran Ice Sheet advanced from the Coast Mountains of British Columbia and covered Vancouver Island and large parts of the current continental shelf (Booth et al., 2003; Clague and James, 2002). The Puget and Juan de Fuca lobes of the Cordilleran Ice Sheet surrounded the Olympic Mountains in the east-southeast and in the north, respectively (Fig. 1b). Alpine glaciers incised deep valleys in the landscape, particularly on the western side of the range (Adams and Ehlers, 2017; Montgomery, 2002), where piedmont glaciers almost reached the Pacific Ocean (Thackray, 2001). Glacial erosion varied across the range, as the location of the Pleistocene equilibrium line altitude increases from 1000 m in the west to 1800 m in the east (Porter, 1964) due to a strong precipitation gradient (> 6000 mm yr⁻¹ in the west, < 1000 mm yr⁻¹ in the east). Determining the exact onset of glaciation in the Olympics has proven difficult, but the oldest deposits of the Cordilleran Ice Sheet in the Puget Lowland are as old as 2 Ma and deeply weathered alpine till on the west side of the Olympics is interpreted to be of the same age (Easterbrook, 1986).

Figure 2Map of new and previously published thermochronometric ages within the Olympic Mountains for (a) AHe, (b) AFT, (c) ZHe, and (d) ZFT. Data from the literature (Batt et al., 2001; Brandon et al., 1998; Brandon and Vance, 1992; Michel et al., 2018; Stewart and Brandon, 2004) are indicated by circles with a black dot. Note that the color coding of the symbols varies between panels. For AFT literature samples, the different reset states (fully reset, incompletely reset, and unreset) are indicated by symbol shape. Maps of exhumation rates as suggested by (e) Michel et al. (2018) and (f) Brandon et al. (1998). The white triangle denotes the location of Mt. Olympus.

Download

2.2 Previous thermochronometry studies in the Olympic Mountains

Within the Olympic Mountains, there is an extensive dataset of thermochronometric cooling ages from bedrock samples (Figs. 1b and 2) for AHe (Batt et al., 2001; Michel et al., 2018), AFT (Brandon et al., 1998), ZHe (Michel et al., 2018), and ZFT (Brandon and Vance, 1992; Stewart and Brandon, 2004). These thermochronometer systems record cooling through a temperature range of ∼60–240 ^∘C (e.g., Brandon et al., 1998; Farley, 2002; Gallagher et al., 1998; Reiners et al., 2004), as they have effective closure temperatures of 70, ∼120, ∼180, and ∼240 ^∘C, respectively, for a cooling rate of ∼10 ^∘C Myr⁻¹ (Ehlers, 2005). The interpretation of thermochronometric cooling ages from sedimentary rocks (such as in the Olympic Mountains) is often complicated when the cooling signal from the sediment source region(s) has not been reset due to reheating during subduction and metamorphism. If a sedimentary rock sample has not had sufficient exposure to temperatures above the closure temperature of a given thermochronometer, the sample might retain cooling ages that represent the source region's cooling history (referred to as unreset) or might be a mixture of provenance cooling histories and the reheating process (incompletely reset sample). Determining whether a sample is completely reset, incompletely reset, or unreset can be difficult and usually depends on the statistics of cooling age populations derived from the dated mineral grains (e.g., Brandon et al., 1998). The reproducibility of single-grain (U–Th) ∕ He ages from a sample provides an indication of whether a sample is reset or not. This is typically determined with n=4–7 grains. For the fission-track method, a larger number of grains is typically dated (n=20–100) to reduce the uncertainty in the final cooling age calculation. For samples with a large population, statistical methods can be applied to decompose the chronometer date distribution into different populations and to determine if some portion of the sample is reset (Brandon, 1992, 1996). In the case in which a sample is incompletely reset, a significant young age peak is determined and interpreted as the sample cooling age (e.g., Brandon et al., 1998).

In the Olympics, the youngest published reset AHe ages (≤2.5 Ma) can be found in the western and central portions of the mountain range, and there are two unreset samples in the east (Fig. 2a). The pattern of AFT ages is more complicated (Fig. 2b), and most reset and incompletely reset samples are located in the central part of the mountain range, whereas unreset samples are restricted to areas outside the central (high-topography) part of the range. The youngest reported AFT samples (2–4 Ma) are incompletely reset samples, and fully reset samples have cooling ages between 7 Ma and 27 Ma. ZHe data show a well-developed trend of unreset cooling ages at the coast and reset 5–6 Ma ages in the headwaters of the Hoh and Elwha rivers (compare Figs. 1b and 2c). Reset ZFT samples (∼13–14 Ma) are confined to a small area east of Mt. Olympus (Fig. 2d).

Based on thermo-kinematic modeling, Michel et al. (2018) attributed the observed AHe and ZHe age pattern to an ellipse-shaped exhumation pattern (with highest exhumation rates in the central, high-topography part of the mountain range; Fig. 2e), as predicted for a mountain range situated in an orogenic syntaxis setting (Bendick and Ehlers, 2014). Here, a bend in the subducted slab creates a mechanical stiffening, which in turn leads to rapid and focused exhumation at the surface (Bendick and Ehlers, 2014). High uplift rates in the central, high-topography part of the mountain range are also corroborated by topographic analyses (Adams and Ehlers, 2017) and denudation rates based on cosmogenic nuclides (Adams and Ehlers, 2018). Furthermore, modeling of particularly young AHe ages (< 2.5 Ma) suggests that exhumation rates increased significantly by 50 %–150 % due to Plio–Pleistocene glacial erosion (Michel et al., 2018).

Table 1Data for the ocean drill cores shown in Fig. 3.

For core ODP 888, information is taken from Westbrook et al. (1994), for ODP 1027 from Su et al. (2000), and for DSDP 174 from Kulm et al. (1973). Sedimentation rates are calculated in this study using the reported thicknesses and age constraints. ^a If total thickness exceeds drilled thickness, then the total thickness was estimated from seismic data (e.g., ODP 888). ^b Due to poor core recovery, the Plio–Pleistocene boundary can only be confined between 418 and 446 m. ^c Ages based on biostratigraphy. For cores ODP 1027 and DSDP 174, the Plio–Pleistocene boundary was recovered and an age of 1.7 Ma is used here as reported by Su et al. (2000). ^d For cores ODP 888 and DSDP 174, the age refers to the age of the oceanic crust and is taken from Fig. 3 at the respective location of the core. For ODP 1027, the age refers to the age of the oldest sediment in the core taken from Su et al. (2000). ^e This rate is calculated for the recovered core interval, which only encompasses 600 kyr. ^f The total sedimentation rate equals the total thickness divided by the age of the oceanic crust.

Download Print Version | Download XLSX

2.3 Offshore sediment record

Data constraining the sediment thickness on the Juan de Fuca Plate before the incorporation of sediment into the accretionary wedge are summarized in Fig. 3. Three boreholes were drilled into the blanketing sediments of the Juan de Fuca Plate during deep-sea drilling projects (ODP 888, ODP 1027, and DSDP 174; Fig. 3 and Table 1) and provide estimates of the sediment thickness and age constraints. The sediment thickness at the deformation front of the subduction zone has been estimated from three seismic studies (Adam et al., 2004; Booth-Rea et al., 2008; Han et al., 2016).

Figure 3Map of the Cascadia Subduction Zone showing the age of the oceanic crust (Wilson, 1993) and sediment thickness, estimated from the sediment cores of ocean drilling programs (holes ODP 888, OPD 1027, and DSDP 174) and seismic studies (Adam et al., 2004; Booth-Rea et al., 2008; Han et al., 2016). The amount of Quaternary sediment material estimated from cores is also included (Kulm et al., 1973; Su et al., 2000; Westbrook et al., 1994); more information about the drill cores is provided in Table 1. The locations of two major submarine fans (Nitinat Fan and Astoria Fan) are indicated by the dotted pattern. The Fraser and Columbia rivers are the main modern sediment sources for the Nitinat and Astoria fans, respectively. White dashed lines indicate the position of cross sections presented in this study (see also Fig. 7).

Download

Most of the sediment is contained within two deep-sea sediment fans with different sediment sources. Today, sediment sources for the Nitinat Fan (offshore of Vancouver Island and the Olympic Mountains) include detritus from Vancouver Island, the Olympic Mountains, and material delivered by the Fraser River system (Fig. 3), which drains large parts of the Canadian Cordillera including the British Columbian Coast Mountains (Carpentier et al., 2014; Kiyokawa and Yokoyama, 2009). The Astoria Fan offshore of the Oregon coast is mostly fed by the Columbia River and is sourced by a large area in the interior of the USA (Fig. 3).

The total sediment thickness varies between 2600 and 3500 m at the deformation front and decreases rapidly to 600 or 900 m approximately 100 km away from the deformation front. At the locations of ODP 1027 and DSDP 174, up to 50 %–70 % of the total sediment thickness is made up of Quaternary deposits, and sedimentation rates more than doubled during the Quaternary (from 80–110 to 250–270 m Myr⁻¹; Table 1). At the location of ODP 888 the drilled 570 m of core was deposited over the past 600 kyr, suggesting very high sedimentation rates of 950 m Myr⁻¹ compared to 400 m Myr⁻¹ for the total sediment thickness of 2600 m at the location of the core (Table 1). As determined from a detailed stratigraphic analysis of core ODP 888, sedimentation rates are also highly variable during the Quaternary. Rates during glacial periods can be as high as 1900 m Myr⁻¹ compared to 700 m Myr⁻¹ during interglacials (Knudson and Hendy, 2009). At sites ODP 888 and 1027, the source region of the sediments has been the Canadian Cordillera for the past 3.5 Myr, which has not been affected by glacial–interglacial cycles (Carpentier et al., 2014; Kiyokawa and Yokoyama, 2009). The provenance of the sediments at DSDP 174 is mostly the Proterozoic Belt Supergroup in the interior of the USA and differs significantly from present-day detritus of the Columbia River (Prytulak et al., 2006). Hence, Prytulak et al. (2006) suggest that the deposition of the upper 630 m of sediment at this site and the buildup of the Astoria Fan were governed by glacial outburst floods.

3.1 Thermochronometric methods

Our strategy with thermochronometric dating was (1) to obtain samples that are multi-dated with up to four thermochronometer systems (because these are particularly sensitive to reveal variations in exhumation rate) and (2) to collect samples within vertical profiles in order to obtain estimates of the exhumation rate at the site of the respective profile. Therefore, we dated several literature samples with additional thermochronometer systems (Table 2) and we also present 19 new bedrock samples from vertical profiles (Fig. 4, Table 2) and two additional bedrock samples (OP1528 and OP1556; Fig. 2, Table 2) collected at an elevation of ∼400 m, enlarging the existing ∼400 m equal-elevation data of Michel et al. (2018). All new samples are sandstones of varying grain size. A sample transect at Mt. Olympus extends from the bottom of the Hoh Valley to the apex of the Olympic Peninsula (Mt. Olympus, 2428 m), covering ∼2 km of relief (Fig. 4a and b). The Mt. Anderson transect starts in the upper reaches of the Quinault Valley and terminates on the flank of Mt. Anderson, covering a total elevation difference of ∼1600 m (Fig. 4a and c). The Blue Mountain transect is located in the northern part of the Olympic Peninsula close to Blue Mountain, covering an elevation difference of ∼1300 m (Fig. 4a and d). All collected samples were dated with the AHe and ZHe techniques; three were dated by AFT, and two were dated by the ZFT technique. Additionally, we dated 13 samples from Michel et al. (2018) by AFT and 5 by ZFT thermochronometry. This process yielded seven samples with AHe, AFT, ZHe, and ZFT cooling ages (Table 2).

Figure 4(a) Map of samples for which three to four different thermochronometer systems are available. The pie charts show the reset stage of a particular thermochronometer system for the sample. If AFT ages are unreset, the peak age of the youngest age population is given as the sample age (see Table 3 for older populations). Ages denoted with an asterisk are taken from Michel et al. (2018). The Hurricane Ridge Fault (HRF) separates the rocks of the accretionary wedge from the surrounding Coast Range Terrane in the hanging wall. The locations of the three different elevation transects (Mt. Olympus, Mt. Anderson, and Blue Mountain) are indicated on the map and the resulting age–elevation plots are shown in (b) to (d). In (b) and (c) the dashed colored lines correspond to possible exhumation rates interpreted from the respective thermochronometer. All uncertainties are 1 standard deviation, SGA: single-grain age.

Download

Standard mineral separation techniques (sieving, magnetic and gravimetric separation) were used to obtain apatite and zircon separates from crushed rock samples. For AHe and ZHe dating mineral grains were handpicked and dated in the thermochronometry lab of the University of Tübingen, following the dating protocol of Stübner et al. (2016). The Ft correction for apatite (Farley, 2002) and zircon (Hourigan et al., 2005) is applied to the measured amount of helium. The (U–Th) ∕ He age equation is solved using the approach of Meesters and Dunai (2005). From each sample, we dated four to seven apatite grains or three to six zircon grains and the results of single-grain analyses can be found in Tables S1 and S2 in the Supplement. Our approach for assessing whether a sample is reset or unreset and the procedure for the exclusion of outliers are explained in the Supplement (Sect. S1.1). For reset samples, we calculate the arithmetic mean age from the accepted single-grain ages, which is reported in Table 2 along with a 1 standard deviation (1 SD) uncertainty.

Table 3Results from fission-track dating.

For AFT and ZFT, 20 grains per sample were dated in a first step and it was checked whether the sample passes the χ² test and can be considered as reset (i.e., >5  %; an indication of belonging to the same age population; e.g., Galbraith, 2005). If so, the pooled ζ age was considered as the sample age and reported here. If a sample failed the χ² test (i.e., <5  %), the sample is considered unreset and, in the case of AFT, 100 grains were dated if enough grains were available. The detrital age distribution was then decomposed into detrital age populations using BINOMFIT (Brandon, 1992, 1996) and the peak ages of those populations (with asymmetric error range for each age peak corresponding to the 68  % confidence interval; CI) are reported. For the ZFT method, information on whether the sample is reset or unreset is sufficient for this study and no further grains were dated. “Fraction” equals the number of grains contained within the respective age peak. N is the number of counted grains. Contrary to Brandon et al. (1998), we did not consider multiply or partially reset AFT samples but treated them as unreset because our thermo-kinematic model can only be applied to fully reset samples. ^a Results for this sample are obtained by merging grains from our sample OP1527 (n=103) and sample AR39 (n=31) from Brandon et al. (1998). ^b Results for this sample are obtained by merging grains from our sample OP1528 (n=68) and sample AR40 (n=12) from Brandon et al. (1998).

Download Print Version | Download XLSX

Fission-track dating of apatite and zircon was performed using an external detector and ζ-calibration techniques (Hurford, 1990). Details about the treatment of the apatite and zircon mounts in the Tübingen thermochronometry laboratory can be found in Falkowski et al. (2014) and Falkowski and Enkelmann (2016). Table 3 contains the AFT and ZFT sample ages and explains the procedure for assessing whether a sample is reset or unreset. Data for single-grain ages from fission-track dating of apatite and zircon are reported in Tables S3 and S4.

3.2 Thermo-kinematic modeling: model setup and boundary conditions

To interpret cooling histories recorded by our thermochronometers as exhumation histories, we used a modified version of the thermo-kinematic model Pecube (Braun, 2003), which contains a built-in Monte Carlo approach to resolve temporal variations in exhumation histories (Adams et al., 2015; Thiede and Ehlers, 2013). The model allows for the exploration of possible exhumation histories for a particular sample by varying exhumation rates through time at defined time steps. The accuracy of a particular exhumation rate history is estimated by comparing modeled with observed cooling ages. More age constraints, and hence thermochronometer systems, lead to better-resolved modeled exhumation histories. Therefore, although we report 21 new thermochronometric ages, we only used the seven samples that have age constraints from AHe, AFT, ZHe, and ZFT in our modeling efforts (OP1513, OP1517, OP1533, OP1539, OP1551, OP1573, OP1582; Table 2).

Table 4List of parameters used for the Pecube modeling.

Download Print Version | Download XLSX

The thermophysical parameters chosen for the modeling are typical values reported for the sandstones of the Olympic Mountains (Table 4). We performed a sensitivity analysis in order to find the most suitable time step for our simulations and the results of that analysis can be found in the Supplement (Sect. S2). Based on the analysis, a time step interval of 1 Myr seems to be most appropriate to use, given the range of our thermochronometry ages and their respective uncertainties. During further modeling, we initiated the models at 20 Ma and used the time step interval of 1 Myr with a maximum testable exhumation rate of 6 km Myr⁻¹. For each sample, we ran 20 000 simulations (each corresponding to a different exhumation history) and assessed the goodness of fit between observed and modeled data for the respective exhumation history using a reduced χ² test. Here, sample ages τ_o were compared with modeled ages τ_m using the uncertainty of the sample age σ_o for the number (N) of thermochronometer systems available for the respective sample:

If χ²≤2, a specific model run was accepted as good. The number of accepted exhumation histories is shown in Fig. 5 for each sample. From the range of acceptable exhumation rates at each time step (shown as blue shaded areas in Fig. 5), we calculated the mean exhumation rate together with 1 standard deviation for each time step (red dashed lines and grey areas in Fig. 5). Although the model provides output for the entire model duration of 20 Myr, a meaningful exhumation rate can only be obtained for the time interval between the oldest thermochronometric age of a sample and today (shown in Fig. 5).

Figure 5Results from thermo-kinematic Monte Carlo modeling for the seven considered samples (OP1513, OP1517, OP1533, OP1539, OP1551, OP1573, OP1582). The location of each sample within the Olympic Peninsula is shown, together with the respective elevation (Elev). The entire range of exhumation rates from the number of accepted model runs (N) is outlined by the blue shaded area, from which the mean rate and 1 standard deviation (1 SD) is calculated at each time step. Black, green, and blue stippled boxes outline measured AHe, AFT, and ZHe ages of the samples with 1 SD.

Download

For our purpose, we focus on exploring temporal variations in exhumation rates and therefore use a 1-D model, whereby each sample is modeled independently from each other. In a 1-D model, heat transport and the movement of particles are only considered in the vertical dimension within a column of rock, ignoring topography. This mode of modeling was selected because it allowed us to efficiently perform thousands of simulations quickly in order to cover a large range of possible exhumation rates. The high number of exhumation histories accurately predicts our observed cooling ages and allows for a robust statistical assessment of the best-fitting exhumation history. Previous publications addressing exhumation histories in other orogens have also highlighted the fact that 1-D models are often sufficient to explain most of the signal recorded in thermochronometric systems (e.g., Adams et al., 2015; Thiede and Ehlers, 2013). In the Olympic Mountains, Michel et al. (2018) argued that exhumation histories for the thermochronometer systems considered here can also be explained well by vertical velocity paths. Because the spatial resolution of our seven considered samples is poor and they are all from the interior part of the mountain range (Fig. 4), we cannot further resolve the exhumation rates outside this area, making a 3-D model very difficult to validate. Therefore, we limit our interpretations to the better-resolved exhumation histories from the 1-D model and focus on the primary temporal changes, rather than paleotopography, or specific differences in the exhumation rates between samples.

Five of the seven considered samples are from the same elevation range (400–580 m), but two samples are from higher elevations (1360 and 1500 m; Fig. 5). Large differences in elevation between the samples can impact the direct comparison between them (e.g., it can affect how changes in exhumation rate are recorded from location to location). However, we are not able to correct for this circumstance (by using an age–elevation relationship) and therefore try to consider this complication when interpreting our exhumation rate histories from the different samples.

Figure 6Constraints used for our quantitative accretionary influx and denudational outflux calculations. (a) Cartoon illustrating our approach for calculating the accretionary influx. The influx corresponds to the sediment scraped off the subducting Juan de Fuca Plate and is governed by the plate convergence rate (u) and the incoming sediment properties (thickness d, porosity η). Because we use a three-dimensional geometry, accretion is considered along a length (l) within a vertical plane. This length corresponds to the length of the coastline indicated in panel (d). OM: Olympic Mountains, CRT: Coast Range Terrane. After Batt and Brandon (2002). (b) Temporal evolution of the plate convergence rate used in the calculations, considering only the component perpendicular to the deformation front (black envelope), and the original output (blue envelope) from the plate reconstruction model of Doubrovine and Tarduno (2008). To provide an uncertainty for our calculations, we consider a range of convergence rates (comprising the width of the envelope) for each time step based on two different rotation models in the model of Doubrovine and Tarduno (2008) (see text for details). (c) Temporal evolution of the plate convergence angle (Doubrovine and Tarduno, 2008) used to correct the plate convergence rate in (b). (d) Exhumation rate pattern from Michel et al. (2018) used for our outflux calculations. The range of displayed rates (0.25–0.9 km Myr⁻¹) corresponds to the rates prior to the glacially induced increase in exhumation rates. The outflux is based on the spatial integration of the exhumation rate pattern. Values for the integrals are listed below the plot for the respective increase in exhumation and time. The white dashed line was used for integrating the exhumation rate using a two-dimensional geometry, which is further explained in the Appendix.

Download

3.3 Methods for estimating flux steady state

To assess the flux steady-state hypothesis of the Olympic Mountains, we need independent estimates of the material influx and outflux over time. For this, we focus on the time period since 14 Ma, which corresponds to the proposed establishment of flux steady state (Batt et al., 2001; Brandon et al., 1998). Flux steady state requires that the material influx into the wedge equates to the amount of accreted material removed from the subducting slab. We assessed the amount of accreted sediment (material influx) with two approaches. First, we calculated the amount of sediment incorporated into the accretionary wedge at the deformation front (Fig. 6a) during the 14 Myr period. Second, we compared this amount of “expected” accreted sediment with the observed amount of sediment residing in the accretionary wedge along two cross sections. The material outflux from the mountain range is estimated using results from thermo-kinematic modeling by equating modeled exhumation with denudation, which can then be integrated spatially and over the 14 Myr period.

Previous flux steady-state analyses in the Olympic Mountains were performed in two dimensions along a profile crossing the Olympic Peninsula. However, exhumation rates within the Olympic Mountains are known to vary spatially (Brandon et al., 1998; Michel et al., 2018). This suggests that the outflux is spatially variable, depending on the location within the mountain range. Hence, we performed our flux analysis in three dimensions and the resulting geometries are summarized in Fig. 6. The influx is calculated along the length of the deformation front, and for the calculation of the outflux we considered almost the entire area of the Olympic Peninsula.

3.3.1 Calculating the accretionary influx

We used a similar approach as Batt et al. (2001) to calculate the accretionary influx, but used a three-dimensional geometry and additionally considered temporal variations in the variables used. Assuming all sediments resting on the subducting oceanic crust are incorporated into the accretionary wedge, the volume of accreted sediment (V_sed) can be approximated using the porosity of the sediment η, incoming sediment thickness d, length of the coast l, the duration of subduction t, and the subduction velocity perpendicular to the present-day deformation front u_per:

$$\begin{array}{}\text{(2)}& V_{\mathrm{sed}}=(\mathrm{1}-\mathit{\eta })\cdot d\cdot l\cdot t\cdot u_{\mathrm{per}}.\end{array}$$

A limitation to this approach is the assumption that all sediment resting on the down-going plate is accreted. There is geochemical evidence that, at early stages of subduction at the Cascadia Subduction Zone, sediment has been incorporated into the mantle and was involved in the magmatism of the Cascades Arc (Leeman et al., 2005; Mullen et al., 2017). However, there are no estimates on the amount of sediment transported into the mantle at present, and most sediments seem to be accreted either at the deformation front or underplated at depth (Calvert et al., 2011).

Table 5Results from influx and outflux calculations using a three-dimensional geometry.

The entire procedure for calculating the influx and outflux is described in Sect. 3.3. The influx volumes are reported as ranges because minimum and maximum convergence rates (Fig. 6b) have been obtained from the plate reconstruction model of Doubrovine and Tarduno (2008). ^a Sensitivity to incoming sediment thickness: the accretionary influx volume is calculated for three different sediment thicknesses, yielding a minimum volume (1.5 km thickness), maximum volume (2.5 km thickness), and a more realistic volume (whereby the volume increases from 1.5 to 2.5 km at 2 Ma). ^b Sensitivity to an increase in exhumation rates: the denudational outflux volume is calculated assuming constant exhumation rates and considering the increase in exhumation rates due to glacial erosion, with an increase by 50 % at 3 Ma or an increase by 150 % at 2 Ma. Exhumation rates are based on Michel et al. (2018) and displayed in Fig. 6d and e.

Download Print Version | Download XLSX

The variable with the greatest uncertainty in this calculation is the sediment thickness back in time that has now been subducted below the Olympic Mountains. As discussed above (Sect. 2.3), the present-day sediment thickness of 2.5 km is the product of increased offshore sedimentation during the Quaternary, and the pre-Quaternary sediment thickness is difficult to determine. Following the approach described in the Supplement (Sect. S3.1), we estimated a pre-Quaternary sediment thickness of 1.5 km. In total, we calculated three different sediment volumes based on different sediment thicknesses (Table 5). Assuming a thickness of 1.5 and 2.5 km for the 14 Myr period yields a minimum and maximum value for the accreted sediment volume, respectively, representing a sediment volume unaffected by Quaternary sedimentation (1.5 km) and a volume for a likely too-high sediment thickness, using the modern thickness (2.5 km). Alternatively, we considered an increase in sediment thickness from 1.5 to 2.5 km at 2 Ma, which likely yields the geologically most meaningful volume.

The porosity of the sediment stack depends on the thickness and decreases with increasing overburden. According to Yuan et al. (1994), the porosity at depth z of the sediment stack can be approximated by

$$\begin{array}{}\text{(3)}& \mathit{\eta }=\mathrm{0.6}\cdot e^{-z}.\end{array}$$

Using this equation, we calculated mean porosities of 31 % and 22 % for our sediment thicknesses of 1.5 and 2.5 km, respectively.

Because the dip direction of the present-day deformation front is 72^∘ (Φ_def) and we only considered accretion perpendicular to the deformation front, we corrected the convergence rate (u) by using the convergence angle (Φ) between the Juan de Fuca and the North American plates:

$$\begin{array}{}\text{(4)}& u_{\mathrm{per}}=u\cdot {\displaystyle \frac{\sin \left(\mathit{\varphi }\right)}{\sin \left({\mathit{\varphi }}_{\mathrm{def}}\right)}}.\end{array}$$

Both the convergence rate and angle are variable over time, and therefore we capitalized on the plate reconstruction model of Doubrovine and Tarduno (2008) to estimate these parameters over the past 14 Myr. Values shown in Fig. 6b and c were calculated using the East–West Antarctica plate circuit model from Doubrovine and Tarduno (2008) for two different rotation models (Farallon M1 and M2 in the original publication). This yields a range of possible convergence rates and angles, providing an uncertainty on the calculated sediment volume. The temporal resolution is given by the number of magnetic isochrons used for the plate circuit reconstruction by Doubrovine and Tarduno (2008). From the temporal evolution of the corrected convergence rate (Fig. 6b), we calculated the sediment volume V_sed accreted during the 14 Myr period using Eq. (2) and the parameters discussed above. For the length l in Eq. (2) we assumed a value of 131 km, which corresponds to the length of the coastline in the area of the exhumation rate pattern (Fig. 6d). The calculated sediment volumes are reported in Table 5.

3.3.2 Sediment volumes along cross sections

We estimated the actual volume of sediment currently residing in the accretionary wedge along two cross sections, which are approximately 50 km apart (Profile 1 and 2 in Fig. 7). The lower boundary of the accretionary wedge is the top of the subducting oceanic plate, which is constrained from the Slab 1.0 model (Hayes et al., 2012; McCrory et al., 2012). The upper boundary is defined by the present-day topography and bathymetry (from 10 and 500 m resolution digital elevation models, respectively) and the Hurricane Ridge Fault (HRF). At the surface, the location of the HRF is adopted from a geologic map (Tabor and Cady, 1978) and below the surface we use information provided by a seismic study at depths of 22 and 34 km (Calvert et al., 2011). The uncertainty related to the position of the HRF (error bars at HRF nodes in Fig. 7) was propagated to estimate an uncertainty for the calculated sediment volumes. A further explanation of this approach is given in the Supplement (Sect. S3.2). Because the location of the HRF is not resolved at greater depths, we truncate the area considered for volume calculation at 34 km of depth. Finally, the calculated volume is corrected for the porosity of the sediment stack. Davis and Hyndman (1989) use porosities of 4 %–10 % for sediments contained within the accretionary wedge offshore of Vancouver Island. Hence, we use an average porosity of 6 % in our correction.

Figure 7Sediment volumes calculated along two cross sections spanning the Olympic Peninsula (Profile 1 and 2, the vertical exaggeration is 2; see inset for location). For an explanation of the procedure used, see the text. The reported uncertainties for the volume are based on the uncertainties in the position of the Hurricane Ridge Fault (indicated with error bars at the respective symbol). Numbers in the inset correspond to the following: (1) position of Profile 1, (2) position of Profile 2, and (3) position of profile by Davis and Hyndman (1989) referred to in the text.

Download

4.1 Thermochronometry

Along the Mt. Olympus elevation transect (Fig. 4b), AHe ages (1.9–3.7 Ma) overlap each other within sample error (except for the uppermost sample). ZHe ages (4.8–8.5 Ma) show a similar behavior (with the exception of the lowermost sample; Fig. 4b). AFT ages for two samples are 5.1 and 6.2 Ma, and the obtained ZFT ages of this transect are all unreset. Within the Mt. Anderson transect (Fig. 4c), AHe ages (1.5–3.9 Ma) increase with elevation up to 1400 m and decrease between 1400 and 2100 m. ZHe ages vary between 6.5 and 8.9 Ma and one sample at ∼1400 m has an AFT age of 7.8 Ma. For the Blue Mountain transect (Fig. 4d), AHe ages (3.6–30.1 Ma, and one unreset sample) do not show a clear correlation with elevation, but, interestingly, the uppermost sample yields the youngest age. ZHe ages of dated samples of this transect are all unreset.

Clear spatial patterns for the multi-dated thermochronometer samples are observable (compare Figs. 2 and 4). AHe ages are reset (apart from one sample in the northeast of the mountain range) and decrease towards the center of the mountain range, where very young ages (< 2.5 Ma) can be found. Seven fully reset AFT samples (5.0–7.8 Ma) are confined to the center of the range (samples OP1513, OP1517, OP1533, OP1539, OP1551, OP1573, OP1582), overlapping the area of reset ZHe samples. The remaining eight AFT samples are unreset (Table 3 and Fig. 4). Two samples at the north and east coast (OP1502 and OP1510) have the youngest age peaks at 26 Ma (comprising 29 % of the dates) and 36 Ma (35 %), respectively. Samples from the western part of the mountain range (OP1521, OP1522, OP1527, OP1528, OP1531) have younger age peaks of 5–16 Ma (comprising 20 %–76 % of the dates). Furthermore, the youngest age peak of these samples decreases in age towards the area of fully reset AFT samples.

We also collected samples (OP1527 and OP1528) close to locations with the youngest AFT ages of Brandon et al. (1998), which were reported as incompletely reset samples (with youngest peak ages of 3.9 and 2.3 Ma). In the original publication, only a small number of grains were dated (n=31 and n=12). To improve the statistics of these two samples, we merge our single-grain ages with those of Brandon et al. (1998) and obtain more robust age distributions (n=134 and n=80; Table 3). The youngest peak ages of the age populations for the two merged samples are 7.4 and 4.7 Ma (2–4 Myr older than age populations reported by Brandon et al., 1998).

ZHe ages constrain an area of reset ages (4.8–10.2 Ma) in the central, high-topography portion of the mountain range (light grey shaded area in Fig. 4a). Five of these samples have AFT (5.1–7.8 Ma) and ZHe (4.8–8.9 Ma) ages that overlap within sample errors, implying rapid cooling (and hence fast exhumation) through both systems' closure isotherms. AHe ages of these samples are younger (1.7–3.9 Ma) and do not overlap AFT ages, indicating that exhumation rates decreased after cooling below the AFT closure isotherm.

Of the seven samples dated with the ZFT method, only sample OP1539 has a fully reset age (12.6 Ma). Together with data from Brandon and Vance (1992) and Stewart and Brandon (2004), this confines reset ZFT samples to a very small area east-southeast of Mt. Olympus, encompassing the headwaters of the Elwha and Quinault rivers (area outlined with a red dashed line in Fig. 4a).

4.2 Exhumation histories from thermo-kinematic modeling

Between 13 000 and 17 800 simulations provide a good fit to the data for each of the seven samples used in the thermo-kinematic modeling (Fig. 5). As expected, the four samples (OP1533, OP1539, OP1551, OP1582; Fig. 5) with overlapping AFT and ZHe ages require fast exhumation rates of > 3 km Myr⁻¹ between 5 and 8 Ma, followed by a reduction to < 0.2 km Myr⁻¹ at 5 or 7 Ma. The reduction of rates for sample OP1573 occurs at ∼9 Ma. However, for this sample the AFT age has a larger uncertainty, and hence we consider the 5–7 Ma decrease in exhumation rates as a more robust signal. Six of the seven samples (except for sample OP1517) also record an increase in exhumation rates at 2–3 Ma to > 1 km Myr⁻¹.

4.3 Estimating the flux steady-state balance

The calculated volumes of the accretionary influx depend strongly on the incoming sediment thickness (Table 5). With our three-dimensional geometry (Fig. 6a) volumes vary among ∼70 000 km³ (1.5 km), ∼76 000 km³ (increase from 1.5 to 2.5 km at 2 Ma), and ∼130 000 km³ (2.5 km). The estimated amount of sediment within the accretionary wedge varies depending on the position within the wedge (Fig. 7). Offshore of Vancouver Island, there is 950–1000 km³ of sediment within the wedge (Davis and Hyndman, 1989), while on the Olympic Peninsula there is up to ∼5300 km³ and 3600 km³ of sediment within the central and southern parts of the mountain range, respectively. Our estimates of the denudational outflux vary for the different exhumation rate scenarios (Table 5) and volumes range from 68 000 km³ for constant exhumation rates to 75 000–82 000 km³ for the exhumation scenario with increasing rates.

5.1 Age–elevation relationships

The cooling ages of samples collected from a quasi-vertical elevation profile (e.g., Fitzgerald et al., 1993; Reiners et al., 2003) can be analyzed by looking at the age–elevation relationship. Often, the purpose is to determine an apparent exhumation rate by fitting a line through the data points when ages are positively correlated with elevation. However, the prerequisite for this approach is that, over the lateral extent of the sampled transect, there is no significant gradient in exhumation rates. This is not necessarily given in the Olympic Mountains (Michel et al., 2018; see also Fig. 2e) and the new data represent this complication (Fig. 4b–d).

At Mt. Olympus, the AHe and ZHe age–elevation relationships do show a positive correlation, suggesting fast exhumation rates of ∼1 km Myr⁻¹ between ∼8 and 2 Ma (Fig. 4b). The Mt. Anderson age–elevation relationship for AHe shows a break in slope at ∼1400 m and decreasing AHe ages at higher elevations, and the large uncertainties of the ZHe ages limit an interpretation (Fig. 4c). While such an “inverse” age–elevation relationship could be caused by a change in relief (Braun, 2002), we interpret it to be a result of the strong spatial variation in exhumation rates along the horizontal distance of the transect (e.g., rates increase from 0.25 to 0.9 km Myr⁻¹ over a horizontal distance of 15–20 km; Fig. 2e). In the case of the Blue Mountain transect (Fig. 4d), we relate the noncorrelation of AHe ages and elevation to an incomplete resetting of the AHe system in this area. Here, some samples experienced high enough temperatures to start, or even complete, resetting of the AHe thermochronometric system, causing the observed variability in AHe ages. All ZHe ages from this transect are unreset, corroborating the fact that this part of the Olympic Mountains has not experienced high temperatures compared to the other transects. Indeed, the Blue Mountain transect belongs to the Coast Range Terrane (CRT), which is at a structurally higher level compared to the accretionary wedge (Fig. 1c). In summary, the age–elevation plots support previous results of strong lateral variations in exhumation and incomplete resetting of thermochronometer systems in the outer part of the mountain range.

5.2 Pattern of exhumation

A well-constrained spatial pattern of exhumation is needed for calculating the denudational outflux. Looking at the spatial distribution of thermochronometric cooling ages provides qualitative information about the pattern of exhumation. In general, the distribution of thermochronometric ages indicates that in the Olympic Mountains the magnitude of exhumation increases from the coast to the center. As discussed above, areas belonging to the Coast Range Terrane (close to the coast or the Blue Mountain area, where unreset AHe ages can be found; Fig. 2a) correspond to the structurally highest parts within the range (Fig. 1c) and were not sufficiently reheated to reset the AHe system. Assuming a geothermal gradient typical for the Cascadia Subduction Zone of ∼20 ^∘C km⁻¹ (Booth-Rea et al., 2008; Hyndman and Wang, 1993) and an AHe closure temperature of ∼60–70 ^∘C, the cumulative exhumation magnitude since the onset of exhumation at ∼18 Ma cannot have been greater than 2–3 km.

The aerial exposure of the accretionary wedge (the Olympic Structural Complex; Fig. 1c) records exhumation from greater depths. Here, all samples yield reset AHe ages, requiring a minimum exhumation depth of 2–3 km. In the center of the mountain range (encompassing the headwaters of the Hoh, Queets, Quinault, and Elwha rivers; Fig. 1b) the area of reset AFT ages approximately overlaps the area of reset ZHe ages (Fig. 4a), requiring deeper exhumation compared to the coastal part of the Olympic Structural Complex.

The area east-southeast of Mt. Olympus (corresponding to the area of reset ZFT samples; Fig. 4a) has been exhumed from the greatest depths within the Olympic Mountains. For an average ZFT closure temperature of ∼240 ^∘C (Ehlers, 2005) and the above geothermal gradient, this corresponds to a maximum exhumation from depths of 10–12 km, confirming previous estimates (Brandon and Calderwood, 1990; Brandon and Vance, 1992).

In summary, the central, high-topography part of the mountain range corresponds to the most deeply exhumed part. This corroborates the exhumation rate pattern (Fig. 2e) suggested by Michel et al. (2018), the pattern of denudation rates based on cosmogenic nuclide dating (Adams and Ehlers, 2018), and results from topographic analysis (Adams and Ehlers, 2017), which all suggest that most of the exhumation–denudation occurs at this location. Hence, we use this pattern for the calculation of the denudational outflux.

Figure 8Summary of volcanic activity, tectonic and climatic events, and convergence rate and angle at the Cascadia Subduction Zone in comparison with our interpreted exhumation rates for the past 25 Myr. Exhumation rates are limited to the time interval covered by our thermochronometric ages (0–11 Ma). The curve depicts the interpreted evolution of exhumation rates based on the modeling results shown in Fig. 5 (see text for details). Volcanic activity after du Bray and John (2011). Tectonic and climatic events are (1) the start of exhumation of the Olympic Mountains (Brandon et al., 1998), (2) the onset of uplift of the Oregon Coast Range (McNeill et al., 2000), (3) rotation in the stress field (Priest, 1990), (4) faster uplift in the Oregon Coast Range (McNeill et al., 2000), (5) Pacific-wide plate reorganization (Wilson, 2002), (6) the onset of North American glaciation (Haug et al., 2005), (7) the onset of glaciation within the Olympic Mountains (Easterbrook, 1986), and (8) a change in the deformational style of the offshore accretionary wedge (Flueh et al., 1998). Convergence rate and angle from Doubrovine and Tarduno (2008).

Download

5.3 Temporal variations in exhumation

Our new thermo-kinematic modeling revealed temporal variations in exhumation rates in the Olympic Mountains (Fig. 5). The decrease in exhumation rates at 5–7 Ma can be readily explained by the reduction in plate convergence rate and the change in convergence direction (Fig. 8). A Pacific-wide reorganization of plate movement at 5.9 Ma has been suggested (Wilson, 2002), and rapid uplift of the Oregon Coast Range at 6–7.5 Ma with a subsequent cessation in uplift has also been attributed to variations in the plate subduction parameters (McNeill et al., 2000). Furthermore, the volcanic record of the Cascadia Subduction Zone shows temporal variations, and the strongest volcanic activity lasted from 25 until 18 Ma (du Bray and John, 2011). A period of volcanic quiescence, lasting from 17 until 8 Ma, was then followed by increased activity starting at ∼7 Ma. A change in the stress field of the Cascadia Subduction Zone occurred at 7 Ma, which likely also affected the composition of the magmatism (Priest, 1990). Therefore, we interpret our observed 5–7 Ma drop in exhumation rates in the Olympic Mountains as a response to changes in the plate tectonic conditions.

In contrast, the increase in exhumation rates at ∼2 Ma indicates a response to climatic rather than tectonic changes. As previously suggested by Michel et al. (2018), increased denudation due to the heavy glaciation of the mountain range led to an increase in exhumation rates by 50 %–150 % starting at 2–3 Ma. Our study corroborates these findings and shows that the observed young AHe ages require a recent increase in exhumation rates from slower rates (< 0.2 km Myr⁻¹) lasting from ∼7 until ∼2 Ma. Glaciation of the North American continent commenced at 2.7 Ma (Haug et al., 2005) and the oldest glacial deposits within the Olympics could be as old as 2 Ma (Easterbrook, 1986), overlapping our modeled increase in rates at ∼2 Ma. Due to the strong spatial variation of the Pleistocene equilibrium line altitude within the Olympic Mountains (Porter, 1964), glacial erosion likely also varied spatially, which could explain the different magnitude of the increase in exhumation rates suggested for the different samples. Increased offshore sedimentation related to glacially eroded sediment affected the deformational style of the offshore wedge, leading to the formation of westward-dipping thrust faults, which changed at ∼1.5 Ma (Adam et al., 2004; Flueh et al., 1998; Gutscher et al., 2001).

Taken together, these observations indicate that temporal variations in exhumation rates within the Olympic Mountains are subject to changes in both the tectonic and climatic conditions (as summarized in Fig. 8). The implication of these variations should be considered for flux steady-state assessment.

5.4 Flux steady state in the Olympic Mountains

5.5 Restrictions and limitations of our approaches

In the sections above, we discussed exhumation in the Olympic Mountains and the results from our flux calculations. In the following section, we want to elaborate on possible restrictions or limitations in our approaches.

With our 1-D modeling, we revealed strong temporal variations in exhumation rates (Fig. 5) related to variations in both tectonic and climatic conditions (Fig. 8). However, two of our modeled samples (OP1513 and OP1517) do not display the decrease in exhumation rates at ∼5–7 Ma. These are from the Elwha Valley (Fig. 4), in contrast to the five samples displaying the decrease, which are located in the western part of the mountain range. This suggests that the response of the orogenic wedge to a variation in the tectonic conditions affects only parts of the wedge and might be controlled by discrete structures. Further sampling and thermochronometric dating would be required to localize possible faults. Furthermore, this places a limitation on the application of a refined 3-D model because it requires constraining parameters such as fault location or displacement on these faults. Besides the importance of single structures, the general pattern of deformation in the Olympic Mountains should still be viewed as controlled by the geometry of the subducted plate (Adams and Ehlers, 2017, 2018; Brandon and Calderwood, 1990; Michel et al., 2018).

Regarding our flux analysis, we based our calculations on the volume of accreted sediment within a certain time (governed by the sediment thickness and the plate convergence rate) and the amount of denuded material (governed by the exhumation rates). As we mentioned in Sect. 3.3.1, a variable with great uncertainty is the sediment thickness over time, which has now been subducted below the Olympic Mountains. In the Supplement (Sect. S3.1) we outline our approach for assessing the pre-Quaternary sediment thickness, which is used in our calculations. Although the reported 1.5 km sediment thickness seems to be a plausible value, we note that this value is afflicted with uncertainties and might have been higher. Nonetheless, our proposed balance between influx and outflux is still tenable if the pre-Quaternary sediment thickness deviated from the assumed 1.5 km. In other words, we suggested an influx volume of 75–78×10³ km³ and calculated outflux volumes between 75×10³ and 82×10³ km³ (Table 5), so even an additional influx volume due to a thicker, unnoticed sediment thickness could be balanced with our calculated outflux volumes. Another simplification in our calculation is the assumption of a spatially uniform sediment thickness over the considered length. Figure 3 shows that the sediment thickness along the deformation front is variable and is highest in the Nitinat and Astoria fans. However, an attempt to reconstruct along-strike variations in sediment thickness over time is challenging and would introduce further uncertainties, and thus we assume an average, constant thickness.

During our influx calculations, we did not distinguish between different modes of accretion, such as frontal accretion or underplating. Batt et al. (2001) concluded that most accretion occurs at the front of the wedge. However, a recent seismic study showed that sedimentary underplating is taking place below the Olympic Mountains (Calvert et al., 2011). For our approach, the mechanism of accretion does not matter because we are only interested in whether mass is balanced over the entire wedge and not at a specific point. As indicated, this is a limitation of our approach and might lead to an overestimation of the actual influx volume because we do not account for the amount of sediment transported towards the mantle.

Flux steady state implies that the outflux from and influx into a mountain range balance each other. An inherent assumption is often that the material removed from a mountain range (the outflux) again enters the mountain range via the influx, which consists of the denuded material from the same source. So in the case of an accretionary wedge, this implies that sediment is recycled and the system behaves as a closed system. As we described in Sect. 2.3 of the paper, the sediment currently entering the accretionary wedge of the Cascadia Subduction Zone is a mixture of sediment from different source regions (e.g., Olympic Mountains, Vancouver Island, Canadian Cordillera, and in the case of the Astoria fan, the interior USA; Fig. 3). With the increased detrital input from the Cordilleran Ice Sheet from outside the Olympic Mountains, this effect became particularly pronounced since the onset of Plio–Pleistocene glaciation. Hence, our influx–outflux calculations for the Olympic Mountains do not represent a closed system, and the influx into the Olympic Mountains is solely controlled by the outflux out of the system. However, our calculations indicate that on long timescales (i.e., over 14 Myr) flux steady state is attained, which might seem surprising given that the sediment thickness is governed by contributions from different source regions. We suspect that processes during sediment deposition, like redistribution by turbidity currents and redeposition in more proximal parts of the Juan de Fuca Plate, play an important role in the final sediment budget. As a consequence, the amount of sediment denuded from the Olympic Peninsula in a given time period (the outflux) is dispersed as it enters the ocean so that for the same time period only a fraction of the sediment thickness (governing the influx) is composed of material originating from the Olympic Peninsula.

Variations in the geometry or extent of the accretionary wedge were also not included in our flux analysis. Since the onset of subduction at the Cascadia Subduction Zone with the present geometry at ∼40 Ma, the wedge must have grown over time in order to attain its present shape. As soon as a balance between accretion and erosion is established, the shape of an orogenic wedge remains constant, controlled by its critical taper (e.g., Davis et al., 1983). However, Adam et al. (2004) showed that the Cascadia accretionary wedge responded to increased offshore sedimentation during the Quaternary through the development of westward-dipping thrust faults, shifting the deformation front further seawards, thereby increasing the extent and volume of the wedge. An important parameter contributing to the shape of the accretionary wedge is the angle of subduction, which is flatter below the Olympic Mountains (compared to areas north or south) due to the bend in the subducted slab (Fig. 1a). A reason hypothesized for bending the subducting slab is extension in the Basin and Range Province starting in the middle Miocene (Brandon and Calderwood, 1990). All these points indicate that parameters controlling the size and volume of the accretionary wedge are both spatially and temporally variable. However, we cannot account for all of these circumstances in our flux calculations because they are difficult to quantitatively constrain from available observations. Furthermore, because we based our flux calculations only on volumes of accreted or eroded material over the 14 Myr period, a comparison between these two volumes should not depend on a change in the shape or extent of the accretionary wedge.

As we pointed out in Sect. 5.4.2, flux steady state is obtained by using a three-dimensional geometry. However, we only considered the deformation front perpendicular velocity component for our influx calculations. The different sediment volumes contained in the reported cross sections (Fig. 7) could indicate that on long timescales additional velocity components must be considered. We can only speculate that margin-parallel transport, which is a contentious topic at the Cascadia Subduction Zone (e.g., Batt et al., 2001; McCrory, 1996; Wang, 1996), also contributes to the accretionary influx. Present-day GPS velocities corroborate this hypothesis, indicating northward movement of coastal areas south of the Olympic Mountains (e.g., McCaffrey et al., 2013; Wells and McCaffrey, 2013).

To summarize, several parameters, like the location of faults within the orogenic wedge, the sediment source region, the temporal evolution of the wedge geometry, and margin-parallel transport, are difficult to constrain from current observations. Although we emphasized that not all of these parameters affect our flux analysis, further knowledge of these will refine the current understanding of steady state in the Olympic Mountains.