Regional exhumation history of brittle crust during subduction initiation, Fiordland, southwest New Zealand, and implications for thermochronologic sampling and analysis strategies

The cooling history of certain minerals can be constrained using thermochronologic analyses that are based upon production and diffusion of noble gases, or rates of fission track creation and annealing. Model cooling ages derived from these observations are commonly used to test thermal predictions from deformation and rock uplift models. While a well-constrained cooling history derived from an individual rock may, in some circumstances, be sufficient to test the hypothesis in question, spatial variations in cooling ages are more usually required to test useful geological models. The spatial distribution of data generally complicates any such analysis because it may not be logistically feasible to collect samples in optimal locations; suitable rock types may not be exposed; relevant data from previous studies may have been collected for different purposes; and limited funds or laboratory facilities usually limit spatial coverage.

The hypotheses we wish to test relate to the process of subduction zone initiation beneath southwestern New Zealand (Fig. 1). The location is globally significant because it is one of the few places on Earth where an active, nascent, and isolated subduction zone can be studied (Gurnis et al., 2004). Oblique convergence during Miocene–Quaternary time has resulted in subduction of Eocene–Miocene ocean crust at the Puysegur Trench and formation of the Alpine fault (Barnes et al., 2005; Christoffel and van der Linden, 1972; Davey and Smith, 1983; Delteil et al., 1996; Lamarche and Lebrun, 2000; Lebrun et al., 2003; Sutherland et al., 2000).

The subducted slab is twisted beneath Fiordland (Fig. 2), with geometry reminiscent of an inverted ploughshare (Christoffel and van der Linden, 1972; Eberhart-Phillips and Reyners, 2001; Reyners et al., 2002). High-amplitude gravity and topographic anomalies exist at both the southern and northern ends of the subduction zone (Fig. 2), and the central part of the zone (47–49°S) is characterized by lower amplitude anomalies and Quaternary subsidence (Collot et al., 1995; Lamarche and Lebrun, 2000). Distributed lithospheric strain that accumulated farther east, within the overriding plate, has bent geologic markers into an orocline, caused crustal thickening and uplift, and resulted in a high-seismic-velocity lithospheric root (Fig. 2) (Eberhart-Phillips and Reyners, 2001; Norris, 1979; Sutherland, 1999). Three young (younger than 5 Ma; probably younger than 2 Ma) and small arc volcanoes exist (Fig. 2), though none are thought to be active (Reay and Parkinson, 1997; Sutherland et al., 2006a).

There is currently a strong spatial correlation between where the subducted slab is being twisted downward and the high topography and positive gravity anomalies of Fiordland (Fig. 2). We hypothesize that this correlation is causally related to the geodynamic setting and that a similar relationship existed in the past (House et al., 2002); hence we use the rock cooling response, caused by exhumation and the emplacement of a slab, and recorded using thermochronology, to track the spatial evolution of the subduction zone.

In our first paper (House et al., 2002) we presented a preliminary data set of 32 apatite (U-Th)/He and 34 apatite fission track ages, showed that thermochonology data could be used to track the evolution of subduction, and presented a hypothesis for how and why exhumation occurred above the developing subduction zone. In our second paper (House et al., 2005) we presented a further 28 apatite (U-Th)/He and 6 zircon (U-Th)/He ages, which increased the geographic coverage, and showed that zircon data provide an essential constraint on the timing of subduction initiation in central and western Fiordland.

In this paper we analyze spatial trends and statistical properties of a much larger data set, which is made up of 66 apatite (U-Th)/He, 210 apatite fission track, 6 zircon (U-Th)/He, and 128 zircon fission track ages (see Supplemental File 1¹); and we also use published K/Ar data. The fourfold increase in the number of data and the inclusion of 128 zircon fission track ages for the first time allows us to reliably determine the subduction-related exhumation distribution and history for the entire region. We consider the statistical properties of the data set and implement a new analysis method: a weighted least-squares regression scheme that allows us to obtain the regional rock uplift history and a quantitative measure of uncertainty surrounding that history. The observations, method, and results have implications for subduction initiation geodynamics, the mechanisms of exhumation and cooling of brittle crust at plate boundaries, and future collection and analysis strategies for thermochronologic data.

The vector v is the direct path between i and j, and the set of path vectors {u_k} is defined by the offsets in equal-age surfaces at each of the k faults that are crossed by v (Fig 3). The term ∇A·ds is the vector dot product between the gradient of the age field and the path vector.

In practice, most sampling is done over a relatively narrow vertical interval, because topographic variations are usually small compared to horizontal distances, i.e., v usually has a shallow inclination. A positive correlation between age and altitude is expected for nearby samples with no faults between them, because isotherms at depth are close to horizontal and ages are frozen into the rock as uplift and exhumation occurs, i.e., ∇A is usually steeply inclined upward; and the magnitude of the spatial age gradient, ||∇A||, is closely related to the cooling rate during the period the age became frozen into the rock, making it a useful quantity to analyze (Benjamin et al., 1987; Ehlers, 2005; Fitzgerald, 1994; House et al., 2002, 1998).

The constant a relates to measurement uncertainty associated with natural variations in helium diffusivity or fission track annealing kinetics (mostly related to fluorine content in apatite), radiation distribution (crystal zoning and inclusions), or morphology of the crystals, or other uncertainty introduced by analytical procedures that is not included in the quoted precision. The term bx_ij is a random walk component associated with faulting and topography. The term c²x_ij² is a first approximation for the age variance introduced by an underlying spatial trend within the age field.

Our data set (Fig. 4; Supplemental File 1 [see footnote 1]) includes 66 apatite (U-Th)/He ages (AHE), 210 apatite fission track ages (AFT), 6 zircon (U-Th)/He ages (ZHE), 128 zircon fission track ages (ZFT), and 8 white mica K/Ar ages from near Milford Sound (Claypool et al., 2002; House et al., 2002, 2005). The K/Ar and U/Pb ages from Fiordland are at least Cretaceous (Nathan et al., 2000), which is older than the time of interest, except in the coastal region near Milford Sound, where published white mica K/Ar ages are Neogene (Claypool et al., 2002).

The age difference and spatial relationships for every age pair of the same type in the Fiordland data set were computed, and results were sorted by separation distance. Based upon previous work (House et al., 2002; Sutherland et al., 2006a; Turnbull and Uruski, 1993), we separated two subgroups of age pairs for each age type on the basis of mean age: ages younger than 10 Ma are inferred to have been frozen into the rock during relatively rapid exhumation and cooling associated with subduction initiation; ages older than 20 Ma are generally inferred to predate the most active tectonic phase and are likely to have cooled at a significantly slower rate through their respective closure temperatures.

Age-altitude correlations (Fig. 5) were assessed by examining age pairs that had relatively small horizontal separations of <3 km and a difference in altitude of >100 m. For samples with horizontal spacing <1.5 km, weak positive correlations are observed for AHE younger than 10 Ma, AFT older than 20 Ma, and ZFT older than 20 Ma. The correlation is noticeably weaker when samples separated by <3 km are considered. The presence of highly variable age-altitude gradients and significant outliers (i.e., age decreases with altitude for some age pairs), combined with the weak (AHE younger than 10 Ma; AFT older than 20 Ma; ZFT older than 20 Ma) or undetectable (AFT younger than 10 Ma) correlation between age and altitude within the data set, suggests that local faulting and block rotations are widespread in Fiordland: the distance between significant faults is inferred to be commonly <3 km.

We note, however, that if we consider the data set as a whole by binning and averaging all data by altitude, then the weak positive correlation observed from individual pairs (Fig. 5: AHE younger than 10 Ma; AFT older than 20 Ma; ZFT older than 20 Ma) can be better quantified (see Supplemental File 1 [see footnote 1]).

A strong linear correlation is observed (Fig. 6) between the square of the age difference and the square of the separation distance (see equations 3 and 4). This is consistent with the observation that the data set displays relatively smooth spatial changes in age over distances of 10–100 km. For AFT and ZFT age pairs with mean age older than 20 Ma, the squared age difference for samples separated by <10 km is significantly underestimated by a linear model, indicating that the bx_ij term (equation 4) is significant at these small offsets. We conclude that b is a function of x_ij: b tends to a constant value at small values of x_ij, but bx_ij tends to a constant value at large values of x_ij (>10 km). We suggest that this is because the true difference in exhumation over large distances is modulated by lithospheric and crustal dynamics and varies smoothly. Note that older ages are from regions that have undergone less total exhumation and are located in eastern Fiordland, farthest from the Alpine fault.

The slope of the linear regression lines is much greater for data pairs with mean age older than 20 Ma than it is for data pairs with mean age younger than 10 Ma (Fig. 6). This follows from previous geological work, which shows that exhumation and cooling rates were much higher during the later period (producing a lower spatial age gradient), and the observation that the slope of the regression line is dependent upon the square of the magnitude of the spatial age gradient, ||∇A|| (see equation 5).

The intercept values of the linear regression lines demonstrate that the geological contribution to age variance between sites is substantially larger than the contribution from analytical uncertainties, even over small distances (<3 km). We can further quantify this “geological noise,” which we define as the geological contribution to age variance, by considering the y axis intercept values and subtracting the contribution from known analytical uncertainty (equations 3–5; Fig. 6). For AHE data pairs with mean age younger than 10 Ma the value is 2.3 My² (Fig. 6), which is equivalent to a mean accumulated throw of 0.5 km, as determined from the mean age gradient of 3.0 My/km (Fig. S1 in Supplemental File 1 [see footnote 1]). For AFT data pairs with mean age older than 20 Ma, a value of 1527 My² and age gradient of 11.8 My/km (Fig. S1) gives a mean throw of 3.3 km. For ZFT data pairs with mean age older than 20 Ma, a value of 770 My² and age gradient of 16.3 My/km (Fig. S1 in Supplemental File 1 [see footnote 1]) gives a mean accumulated throw of 1.7 km. This mean accumulated throw is an estimate of the standard deviation of the sum of the throws on all the faults between two sites, after the throw associated with continuous variation is removed. Although these independent calculations yield similar results, we reiterate that the uncertainties are large, because correlations between age and elevation are highly uncertain (Fig. 5) and dependent upon a number of assumptions that are dealt with more thoroughly in the following.

By taking age pairs with horizontal separation <10 km, we are able to analyze the variance of age pairs as a function of mean age. We have plotted the square root of the mean squared age difference (RMS) as a function of age (Fig. 7), because we infer (from equation 5), that this quantity is linearly related to the magnitude of the spatial age gradient, and hence approximately inversely related to the exhumation rate at the age of cooling. The AFT data set is large and many AFT ages span the time of interest, so it was subdivided into north, central, and southern groups.

ZFT data have high RMS age difference for times before 18–20 Ma (Fig. 7), suggesting slow rates of cooling before 20 Ma. ZFT ages of ca. 16 Ma are most prevalent in the Doubtful Sound region. A small group of ZFT samples from near Milford Sound have elevated RMS age difference values with mean age 12–14 Ma, suggesting a later onset of rapid cooling in that region. The time of onset of low AFT RMS age difference values becomes younger toward the north, suggesting a progressive northward onset of rapid exhumation over the period 16–6 Ma.

If exhumation rate were spatially variable, but not temporally variable, and geological noise were similar everywhere, then the RMS age difference values should plot approximately on a straight line that passes through the origin. It is clear from Figure 7 that many apatite data are close to a straight line that passes through the origin; the slope of the line is 0.3–0.5. Assuming a closure depth of 3–6 km, this implies an accumulated mean throw between samples in the range 0.9–3.0 km. However, it is also clear that all of the data cannot be solely explained by geological noise and temporally-constant, but by spatially-variable rock uplift.

It is clear from our spatial analysis of the data set (above) that a coherent spatial signal over distances of >10 km is embedded in the data set, and that analytical uncertainty is not sufficient to explain observed age variations over even quite small separations of <3 km. We wish to find a model of regional rock uplift history that best reproduces a set of n cooling ages {A_i} of different types and known standard errors {σ_i}, and we would like an estimate for the range of possible models. We choose to implement a weighted least-squares regression scheme that takes into account both analytical uncertainties and geological noise; the model parameterization and assumptions are tailored to the geological hypotheses being tested and the data distribution.

The scheme proceeds as follows (Fig. 8): a rock uplift history is defined from a set of parameters (initially guessed) at each node within a regular spatial grid; a rock uplift history is interpolated at each datum; a thermal history is computed at each datum; a model age is computed at each datum; a residual rock uplift value and effective cooling depth are determined at each datum by perturbing the rock uplift history (and hence thermal history) until the model age matches the observed age (Fig. 9); a parameter search is conducted to find a new regional uplift history that better matches the effective cooling depths and cooling ages; the original parameters are adjusted to be closer to the best parameters found during the search; and the procedure is restarted. The scheme is complete when no improvement (at some specified precision) is found during the parameter search.

We choose a parameterization for the rock uplift history that suits the hypotheses we wish to test: how much total rock uplift has occurred; when did uplift start; and how has the rate of rock uplift changed during the past 5 My? We require three parameters: an initial depth, an age of inception (older than 5 Ma) for rock uplift and exhumation, and a rock uplift value at 5 Ma. A continuous rock uplift history is produced by linear interpolation and then smoothed by recomputing the gradient over a 5 My window of width (Fig. 9). Each rock uplift history starts at 80 Ma and is constrained not to allow burial (subsidence). Full details of the regression scheme and inference of confidence limits are given in  Appendix 1.

The scheme that we choose to implement is purposefully naive with regard to boundary conditions that we impose. This is because our primary objective is to identify the general pattern of rock uplift that our data demand, rather than impose a spatially and temporally variable set of boundary conditions that we expect will correlate with the exhumation and cooling signal that we can measure. The sensitivity of our model and hence conclusions to more realistic boundary conditions are discussed after presentation of our results.

The maximum rock uplift during the time of interest (younger than 34 Ma) was in northern Fiordland, where ~15 km of uplift is inferred (Fig. 10); peak temperatures were higher than the zircon annealing temperature of ~260 °C (Tagami, 2005), and close to the K/Ar white mica closure temperature of ~350 °C (Dodson, 1973; Hames and Bowring, 1994). Farther southwest in coastal Fiordland, total rock uplift reduces to ~12 km between Doubtful Sound and Dusky Sound (Fig. 10), where a zircon partial annealing zone meets the coast. The southeastern margin of Fiordland has undergone only ~4–5 km of uplift during the time of interest (Fig. 10), approximately corresponding to exhumation from the annealing depth for apatite fission tracks.

The onset of rapid uplift was later toward the northeast of the northwestern (Alpine fault) margin of Fiordland (Fig. 10). Uplift started near Dusky sound during the interval 25–15 Ma; reached a maximum near Doubtful Sound during the interval 15–10 Ma; continued to migrate past George Sound during the interval 10–5 Ma; and has been most intense near Milford Sound during the interval 5–0 Ma (Fig. 11; Animations 1 and 2). As the zone of maximum uplift rate spread northeast, a wave of increased uplift rate also spread eastward across Fiordland during the interval 10–0 Ma (Fig. 11). Optimal models for southwestern nodes at Dusky Sound and Doubtful Sound predict that uplift rates have reduced during the past 5 My, but the difference in rate is barely significant, or is not significant (Fig. 10).

Model ages are predicted for every point on the existing topographic surface (Fig. 12), and the inferred thermal structure of the shallow crust (0–1 km below sea level) has a zone of elevated geothermal gradients where Pliocene–Quaternary uplift rates are predicted to be highest (Fig. 13).

The confidence limits shown in Figure 11 are based upon chi-squared tests and hence represent an estimate of the real range of rock uplift histories, because ϵ_ij² >> σ_i² + σ_j² (equation 3; Fig. 6). In general, a model of tectonic stability followed by a phase of high uplift rate (~1 mm/yr) adequately fits most data (Fig. 10). However, spatial mismatch between the smooth model and geological reality is confirmed by examination of rock uplift residuals.

A systematic spatial correlation between rock uplift residuals is most clear for ZFT samples (Fig. 14): the optimum regional model of rock uplift systematically underestimates the true uplift for samples in an ~20-km-wide strip extending from the head of Dusky Sound to Milford Sound. A systematic spatial correlation is also observed for AFT samples in a 30-km-wide region south of Milford Sound, and in a region of 20 km width in the inner parts of Dusky Sound (Fig. 14). The 25:50 km node spacing and curvature limitation imposed during spatial interpolation make it impossible for the smooth model to fit these data.

Where AFT and ZFT residuals of opposite sign occur (inner Dusky Sound and just south of Milford Sound; Fig. 14), then the model limitation of having only one possible inflection point at 5 Ma may also contribute to systematic misfit. A local region of rapid uplift at, for example, 20–10 Ma followed by a slower uplift rate cannot be fit using our chosen parameterization. It could also be that systematic inaccuracies associated with annealing or helium diffusion kinetics (e.g., related to rock type) could contribute to temporal misfit, even given a model with constant uplift rate.

The depth residuals from data that were used to constrain the rock uplift history are shown in Figure 15. Each residual is the vertical perturbation from the best-fitting regional rock uplift model that exactly reproduces the measured age at a sample site (Fig. 9). The combined set of residuals fits well a normal distribution with a mean of 0 km and standard deviation of 1.3 km (mean s_g value). The most striking feature of the distribution of residuals is that there is a systematic bias to positive residual values for systems with high closure temperature, and a slight bias toward negative residuals for systems with lower closure temperatures (Fig. 15). This bias is an inevitable consequence of the fact that low-temperature thermochronometers detect anomalously low exhumation histories, and high-temperature thermochronometers are biased toward anomalously high exhumation histories (see Fig. 16). We acknowledge that this introduces statistical bias into our modeling procedure, but reiterate that the combined residuals are approximately normally distributed.

The influence of faults is resolvable at a local scale in only a few instances. The clearest example is at the southern node (Lake Hauroko), where data clearly indicate a divergence in uplift trajectories: one path records rapid exhumation starting ca. 8–10 Ma, while the other path suggests a later onset ca. 2–4 Ma (Fig. 10, node 0). Examination of sample locations in the context of a geological map reveals that the variable signal is associated with exhumation on either side of the significant mapped Hauroko fault (Turnbull and Uruski, 1995).

Field mapping undertaken during this project and in other regional studies indicates that many brittle faults are recognizable in Fiordland, but no description of the regional pattern of brittle faulting has been published; many faults have only small throw, are difficult to recognize or trace over any significant distance, and kinematic indicators or offset markers are rare. Assuming a hypothesis that many brittle faults exist, our statistical consideration (section 3; Figs. 6 and 7) suggests that the accumulated throw on faults between sites spaced <10 km is in the range 0.5–3.3 km. The results of our bestfit model (Fig. 10) indicate consistent geological noise values in the range 0.6–2.3 km. These values represent accumulated throw on several faults and include limitations imposed by model curvature, so individual fault throws are most likely to be lower than the values given above. Geotechnical data acquired during construction of the Manapouri tunnel, which is used to transport water from Lake Manapouri to Doubtful Sound (Fig. 4), confirm that the rock mass is pervasively fractured and faulted, with many recognizable faults at <1 km spacing.

We adopt an objective but simplified approach to convert rock uplift histories into thermal histories and hence model ages. Specific issues that require discussion are: initial conditions; heat productivity in the crust; the upper thermal boundary condition; edge effects caused by trench topography; and deformation of the lithosphere, particularly the effect of a subducted slab on the basal temperature boundary condition.

We assume a constant initial geothermal gradient that is consistent with that measured in nearby petroleum wells that are outside the zone of intense Cenozoic deformation (Funnell and Allis, 1996). Initial formation of the plate boundary through New Zealand in Eocene–Oligocene time was associated with minor rifting near Fiordland (Norris and Turnbull, 1993; Sutherland, 1995; Turnbull and Uruski, 1993), and may have perturbed the thermal structure of the lithosphere. The effect of locally raised heat flow values may cause our method to overestimate the total exhumation in some places, but we have no basis for imposing initial conditions that are other than uniform.

We ignore heat productivity within the crust, which in Fiordland is mainly composed of dioritic amphibolite or granulite gneisses that were exhumed from the lower crust during Cretaceous time (Clarke et al., 2000; Gibson and Ireland, 1995; Mattinson et al., 1986; Oliver, 1980). Based upon geochemical considerations, Fiordland crust is similar to or less radiogenic than that found in Taranaki Basin (Fig. 1), and hence contributes <18 mW m⁻² to the surface heat flow (Armstrong and Chapman, 1998). Heat productivity within the upper 5–15 km of crust, which is the part of interest, probably accounts for <10% of the background ~60 mW m⁻² of heat flow. We hence conclude that variations in crustal heat productivity are unlikely to alter the geothermal gradient in the upper crust by more than ~1–2 °C km⁻¹; the perturbation caused by rock advection, which we model, is 5–20 °C km⁻¹ (Fig. 13). In addition, the effect is spatially uniform and introduces a comparable or smaller uncertainty than that of the initial conditions. Although this may affect absolute values of rock uplift, our conclusions are not likely to be affected because they are based upon the spatial pattern of rock uplift histories and are not sensitive to small differences in the total exhumation.

We assume a constant upper temperature boundary condition at sea level that is close to the mean annual temperature. Our justification comes from rock and groundwater inflow temperature measurements made during construction of the Manapouri tunnel (Fig. 17), which show our assumption to be a good approximation in central Fiordland.

It is clear from Manapouri tunnel temperature data that groundwater efficiently removes heat from topographic anomalies, which have typical wavelengths of <10 km in Fiordland (Fig. 4), and maintains an almost constant rock temperature at sea level. We observe a greater effect in the western end of the tunnel, which is proximal to the Alpine fault and has greater total Cenozoic exhumation. It appears from the temperature data that the western end of the tunnel is hydrologically distinct from the eastern end, but both are strongly affected by groundwater heat advection (Fig. 17). We are undertaking an additional three-dimensional (3D) modeling study to fully evaluate the implications of data from the Manapouri tunnel.

In the European Alps, observations suggest that significant advection of heat by ground-water is unusual and limited to regions around confined fractured aquifers (Kohl et al., 2001; Marechal et al., 1999; Pastorelli et al., 2001). Hence, we suggest that Fiordland defines one end-member type of mountain range, because the rock mass is extensively fractured (above), there is no sediment or ice to seal the topographic surface, and the region has very high rainfall, which varies from 3.6 m/yr at Lake Manapouri to >8 m/yr on the mountain tops.

There is a topographic slope that falls ~4 km over a distance of ~20 km northwest of Fiordland, where the most western thrust faults deform and accrete basinal sediments (Figs. 1 and 5) (Barnes et al., 2002). The topographic slope and accretion of sediment along the northwest margin of Fiordland require lateral heat flow and will cause our model to overestimate temperatures at depth in the west. Present-day heat flow offshore along the northwest margin of Fiordland is ~40 mW m⁻², as inferred from bottom-simulating reflectors (Townend, 1997), significantly lower than the value of ~60 mW m⁻² determined from four petroleum exploration boreholes just east of Fiordland (Funnell and Allis, 1996).

We accept that the constant temperature boundary condition we assume at 30 km depth is naive, and is probably the most significant issue that makes our results unrealistic. We choose a simple and clearly stated basal boundary condition because it is objective and does not introduce any regional signal; it is our primary goal to recover the regional signal, so we do not want introduce circular reasoning by imposing a particular regional mechanical model.

There are two main issues that affect the realism of the lower boundary condition for our thermal model: emplacement of a slab in the mantle, and deformation within the crust. Whatever the reality, the thermal effects are insignificant to start, but become increasingly important as the deformation process continues. Fortunately, we are most interested in the start of the process, i.e., subduction initiation, so we contend that the time and place of initial cooling and the original depth are reliably determined (although see comments above).

Consider the effects of slab or crustal root emplacement on a particle that is initially in the middle crust at 5–15 km depth, the depth range in which we are interested. Both scenarios advect heat downward and lead to the lowering of temperature gradients and temperatures in the middle crust. This cooling effect from below is in the same sense as that caused by the exhumation of rock from above. Hence, if our premise holds that both effects start at about the same time, then our model rock trajectories are likely to overestimate the rock uplift rates that immediately follow subduction initiation, and must correspondingly underestimate the most recent rock uplift. We suggest that this, in combination with edge effects, is a significant cause of the underprediction of Quaternary uplift rate at the western edge of Fiordland (Fig. 11), where uplifted marine terraces show the rate to be 0.5 mm/yr (Kim and Sutherland, 2004; Turnbull et al., 2007).

The northeast migration of rapid exhumation and broadening of the exhumed zone have previously been suggested from consideration of a subset of the same data (House et al., 2002, 2005). We are now able to quantitatively track the development of this zone using the entire data set (Fig. 11; Animations 1 and 2) and estimate the total amount of exhumation. All directions referred to below are in the Pacific plate frame of reference.

The zone of maximum exhumation rate migrated northeast along the Alpine fault margin of Fiordland at ~1 cm/yr during the interval 20–5 Ma. This rate is ~30% of the relative plate motion rate (Fig. 18) and <50% of the Quaternary Alpine fault slip rate (Sutherland et al., 2006b). Therefore, we conclude that the signal is not caused by an initial condition associated with the geometry of the Australian plate and/or an asperity of the Alpine fault (Barnes et al., 2005; Sutherland et al., 2000). Instead, we suggest that Fiordland exhumation, inferred to be associated with elevated topography, was created in dynamic response to the twisting downward of the obliquely subducted slab. We propose that this zone of downward bending migrated northeast and spread out to the southeast as the overriding plate was internally shortened and tectonically eroded, and the slab increased in areal extent. The tectonic shortening and erosion of the overriding lithosphere have resulted in a lithospheric root with high seismic velocity (Eberhart-Phillips and Reyners, 2001), oroclinally bent basement geology (Norris, 1979; Sutherland, 1999), and a broad (250 km) zone of crustal thickening and mountain building to the north and east of Fiordland (Fig. 18). The coincidence of high-amplitude gravity and topographic anomalies with the most tightly folded part of the subducted slab and the current zone of highest uplift rates in Fiordland (Fig. 11) lends support to our hypothesis that the topography and hence exhumation anomaly is primarily associated with inception and downward deflection of a slab.

As plate motion changed and total convergence increased through time, the subduction interface grew in extent. The zone of downward twisting migrated northeast and the southern tip-line of the subduction thrust migrated southwest. Stress levels are inferred to have been greatest at both these ends, consistent with the current locations of the highest-amplitude gravity anomalies (Fig. 2). It may also be that subduction of fracture zone weaknesses has influenced the geometry of the slab and overriding plate (Lebrun et al., 1998; Malservisi et al., 2003), or that membrane stresses within the slab have an additional effect (Reyners et al., 2002). The zone of subsidence that exists in the wake of the southward propagating subduction thrust tip may be due to some combination of the factors: tectonic erosion of the overriding plate; slab pull; loading of the overriding plate by elevated topography at the northern and southern ends of the zone; changes in bulk-rock or fault-rock physical properties through time, as total displacement increases; and an increase in age of the subducted slab and overriding plate through time—the subducting oceanic lithosphere is currently 40–20 Ma, but local convergence has occurred since 25–15 Ma. The last mechanism is an inevitable direct consequence of initiating subduction at a young fracture zone.

Dynamic support of ~1–2 km of anomalous surface topography is predicted by numerical models to balance the integrated shear stresses associated with subduction initiation processes at depth (Toth and Gurnis, 1998). We acknowledge that more sophisticated 3D models will be required, but predict that a similar result will be found. We note that previous 3D models of Fiordland have suggested that a tear in the subducted slab may be necessary to explain the observed topography and gravity (Malservisi et al., 2003), but we suggest, based on observations from the Kermadec Trench (Billen and Gurnis, 2005), that a viable alternate possibility is that the elastic strength of the subducted lithosphere is much lower than that of intact lithosphere at the surface.

Assuming that a dynamic response of 1–2 km is physically reasonable, then it is expected that material will be removed by erosion from the elevated upper plate, and that an isostatic response to exhumation will restore the topography. In the absence of crustal strain this will cause the Moho to be advected upward. Using this reasoning, we previously calculated that the existing gravity and topography are consistent with ~7 km of crustal thinning, at which point an isostatic balance is reached. In central Fiordland (nodes 4–7), our data are broadly consistent with this hypothesis. However, we estimate 10–15 km of exhumation in western Fiordland (nodes 8–11).

The most likely explanation for the greater value of exhumation than expected from dynamic uplift alone is that crustal shortening has also occurred. Assuming an approximate width of the anomalous zone of 40 km, anomalous exhumation of 6 km, and initial crustal thickness of 30 km, a minimum shortening of ~8 km is calculated assuming conservation of crustal volume. This minimum value may significantly underestimate the true value of shortening if a conveyor of material is coming though the system along a dipping structure, as has been the case farther north along the Alpine fault (Walcott, 1998). However, the Alpine fault is a steeply dipping strike-slip fault with little or no Quaternary convergence across it adjacent to Fiordland (Barnes et al., 2005). Therefore, we conclude that the cause of the exhumation is primarily dynamic uplift associated with subduction initiation, but that there is a significant component of crustal shortening, especially in the west. We cannot rule out the possibility that the Alpine fault may have accommodated significant shortening adjacent to Fiordland during its Miocene phase of activity.

A large body of published literature concludes that in many places the thermal structure of the upper crust is primarily controlled by conduction of heat, and hence low-temperature thermochronologic data can be used to directly derive information about the growth of finite-amplitude topographic anomalies (Braun, 2002; Kohl et al., 2001; Stuwe et al., 1994). However, data from the Manapouri tunnel (Fig. 17) demonstrate that there is no correlation between elevated sea-level temperatures and overlying mountainous topography in Fiordland. In fact, a negative correlation is found if data from the same hydrologic zone are considered. Hence, we conclude that it is not straightforward to determine any information about the evolution of individual topographic features in Fiordland using thermochronologic data. We have started a separate 3D modeling study to investigate what can be learned from the Manapouri tunnel data about the permeability of the crust, and further work is required to assess the viability of learning about topography development from thermochonologic studies in fractured and wet environments like that of Fiordland.

It may be that the onset of glaciation since 2 Ma has had some regional effect on our model by both cooling the rock mass (surface temperature boundary condition) and through enhanced erosion and hence exhumation. However, all of Fiordland has been heavily glaciated, so we might expect any effect to be uniform across the region and not affect the regional patterns in which we are most interested. In addition, most of the cooling ages that we model predate the onset of glaciation, so we expect the effect to be relatively minor. Our model is not designed to test hypotheses concerning changes in exhumation rate due to the onset of glaciation.

We note that there is a correlation and causal link in South Island between the long-wavelength (>50 km) topographic elevation, the local amplitude of short-wavelength (<20 km) topography, and the local exhumation rate. The strength of the correlation is modulated by regional factors such as the rock mass strength, precipitation, and proximity to the coast.

Age-elevation relationships have been widely used to constrain past exhumation rates (Benjamin et al., 1987; Fitzgerald, 1994; House et al., 2002, 1998). The principal assumptions required are that equal-cooling-age surfaces are close to horizontal now and when they were frozen into the rock, and that there are no faults between the sample sites. These may be valid assumptions or a useful approximation in many cases, especially where exhumation rates are low and uplift is driven by processes that do not involve widespread fault movement. However, we suggest that these assumptions do not hold true in an actively deforming obliquely convergent orogen with high upper crustal strain rates, as we find in Fiordland.

Age-elevation thermochronological transects require moderate or high topographic slopes, so that a range of elevations can be sampled as nearby as possible. However, it is well known that active faults produce fault scarps that may grow into mountains, and hence we expect fault locations to be spatially correlated with high topographic slopes and high elevations. Brittle fault zones, especially immature zones with relatively low total displacement (<2 km), are typically complex and composed of numerous strands with rotated blocks between them, or they are associated with zones of distributed folding. It is a significant problem that individual brittle faults are typically manifest at the surface as very narrow zones (<20 cm) that are preferentially eroded and covered by sediment and vegetation, and are hence difficult to identify from imperfect exposure in the field. We hypothesize that the steep locations that superficially seem best for age-elevation transects are affected in many places by brittle faults, and that it can be very difficult to demonstrate the existence of (or lack of) faults from field observations.

We did not design our sample strategy to resolve the details of local faulting, but we sampled several age-elevation transects, and attempted to collect nearby samples with maximum elevation difference in locations where it was practical to do so. When we compare age determinations from nearby (<1.5 km, <3.0 km), we find a poor positive correlation between age and elevation (Fig. 5). We suggest that this results from spatial variability between two competing effects: (1) erosion in the absence of local faults, and (2) active growth of reverse fault zones that lift rock more rapidly in their upthrown mountainous hanging wall. Effect 1 produces the conventional positive correlation between age and elevation, and effect 2 may produce a negative correlation if the difference in exhumation rate is sufficiently large across the reverse fault zone. Our proposal of numerous brittle faults in Fiordland with relatively small offset (<2 km) is consistent with field observations made by us and our colleagues, and with a detailed geotechnical study of the Manapouri tunnel; but we acknowledge that our data set is too sparse over most of the region to demonstrate the existence of any specific faults.

Our statistical analysis of age difference as a function of horizontal separation (Fig. 6) shows that within each thermochronological system there is a coherent regional horizontal age gradient that is likely controlled by lithospheric processes associated with the plate boundary. It is also clear from Figure 6 that the y axis intercept is resolved in the sense that differences in age between nearby samples are much greater than can be explained by analytical uncertainty. However, the intercept values (Fig. 6) do not distinguish between different reasons for local age variation (effects 1 and 2), and amalgamate data from many localities using an unrealistic (planar) function for the local geometry of age variation.

Our regression model provides an alternative measure of local short-wavelength brittle faulting around each grid node, because the regional exhumation model removes long-wavelength lithospheric-scale exhumation signals. The parameter s_g (equation 8) is a summary statistic (standard deviation, km) that quantifies the misfit in rock uplift computed for localities closest to each node, and takes into account the age uncertainty. We suggest that offsets on local faults are a significant contributor, but acknowledge that spatial correlations between residuals (Fig. 14) imply some of the misfit is related to intermediate-wavelength signals that cannot be described by the small number of parameters we choose for the regional model.

Our regional model fits the data using conventional reasoning of subhorizontal isotherms at depth and uniform uplift, and predicts a positive age-elevation relationship everywhere. The result is that the rock uplift residuals have had the regional tectonic signal and the conventionally predicted age-elevation relationship removed. If the true rock uplift on mountains is higher than in valleys, because of the reverse faults that we hypothesize, then we would expect a positive correlation between rock uplift residuals and the elevation of the sample. In Figure 18 we analyze the rock uplift residuals from the regression, by binning them by elevation range of the sample, and then taking the mean value and standard error of each bin. The figure shows that there is a positive correlation, on average, when the entire data set is considered. However, the exact statistical meaning of Figure 18 is difficult to assess because we have amalgamated data of variable spatial density and closure temperature that have already been fit by a regression method. Our data set is barely dense enough to address this issue, and we do not have the computational resources to (marginally) improve our statistical treatment (requires running the regressions many times to make predictions independent of the sample). We certainly cannot resolve details of individual faults around specific topographic features. However, on the basis of Figure 18, we conclude that there is a significant signal embedded in the residuals that is consistent with our hypothesis, and it could be investigated and quantified further.

Our results enable us to track the progress of subduction initiation since 25 Ma in Fiordland, New Zealand. This is significant because there are few other places on Earth where it is possible to observe and analyze this process in such detail (Gurnis et al., 2004). We observe the onset of rapid exhumation at 25–15 Ma in southwest Fiordland, immediately following a time of significant change in regional plate motions (Cande and Stock, 2004). During the period 15–5 Ma, the locus of rapid exhumation broadened and migrated toward the northeast at ~30% of the plate motion rate, but exhumation remained localized along the northwest margin (Fig. 11). Since 5 Ma, the zone of rapid exhumation has become broader (expanded toward the southeast) and the present high-amplitude gravity and topographic anomalies are spatially associated with the most tightly folded part of the subducted slab.

We infer that the exhumation history tracks the development of topography, which is in turn related to the development of a through-going subduction interface and deflection of the nascent slab downward. The zone of subsidence southwest of Fiordland is in the wake of the southward-propagating subduction thrust tip. Subsidence may be due to a combination of factors such as tectonic erosion of the overriding plate, slab pull, loading of the overriding plate by elevated topography at the northern and southern ends of the zone, changes in bulk-rock or fault-rock physical properties through time, as total displacement increases, and an increase in age of the subducted slab and overriding plate through time. Oceanic lithosphere being subducted now has an age range of ca. 40–20 Ma, so it was very young during subduction initiation at 20 Ma.

We conclude that as plate motion changed and total convergence increased through time, the subduction interface grew in extent, a lithospheric root with high seismic velocity was formed to the northeast of the slab (Eberhart-Phillips and Reyners, 2001), and the crust of the overriding plate was shortened and bent into an orocline (Norris, 1979; Sutherland, 1999) (Fig. 19). The zone of downward twisting progressively migrated northeast and the southern tip-line of the subduction thrust migrated southwest. Stress levels are inferred to have been greatest at both these ends, consistent with the current locations of the highest-amplitude gravity anomalies (Fig. 2). Subduction of fracture zone weaknesses may have influenced the geometry of the slab and overriding plate (Lebrun et al., 1998; Malservisi et al., 2003), or membrane stresses within the slab may have an additional effect (Reyners et al., 2002).

We find local positive correlations between age and elevation within our data set, and we find locations where the opposite is true. Data from the Manapouri tunnel show little variation in sea-level temperatures in highly fractured rock beneath high (>1 km) topography. We suggest that our observations are best explained by brittle faults, and conclude that faulting can be significant during regional exhumation at temperatures <300 °C at active plate boundaries. We suggest that an entire population of brittle fault throws and block rotations is likely to exist in such settings, even if their existence is hard to prove. Faults may displace and rotate equal-cooling-age surfaces, and fractures increase rock permeability and encourage thermal perturbation of the upper crust by groundwater flow.

Our study highlights the importance of direct determination of past cooling rate on the same sample, or on very closely spaced samples with known field relations (outcrop rather than mountain scale), rather than relying on age-elevation relationships. A spatially complete and uniform sampling strategy is required to gain insight into the regional pattern of exhumation, which is controlled by lithospheric-scale processes, and the local controls of faulting on topography, where we expect causal relationships and spatial correlations between topographic features and fault locations. Our regression method, although naive in its assumptions, allows us to simultaneously consider the implications of all available data, evaluate summary statistics that describe the geological noise, and interpolate our regional exhumation model to make cooling age predictions that can be tested at any point. The residuals from our regression indicate that, on average (Fig. 18), mountain tops in Fiordland have undergone slightly greater rock uplift than adjacent valleys, even though we cannot identify specific faults.

We thank Lawrence Gaylor, Xu Ganqing, Ian Whitehouse, and Ivan Liddell for fission track sample collection and technical assistance in the University of Waikato Fission Track Laboratory; K. Gallacher and T. Ehlers for providing software that calculated model thermochronologic ages from thermal histories; Simon Cox and Andy Nicol for reviews of an early draft; and Randy Keller, David Shuster, and two anonymous reviewers. This paper is contribution 10,011 of the Division of Geological and Planetary Sciences and 94 of the Tectonics Observatory, Caltech. It was funded by the NZ Marsden Fund and the US National Science Foundation.

The rock uplift since geological age t at grid node g is R_g(t). We choose a parameterization for the rock uplift history that suits the hypotheses we wish to test. How much rock total rock uplift has occurred; when did uplift start; and how has the rate of rock uplift changed during the past 5 My? We require three parameters: an initial depth, an age of inception (older than 5 Ma) for rock uplift and exhumation, and a rock uplift value at 5 Ma. A continuous rock uplift history is produced by linear interpolation and then smoothed by recomputing the gradient over a window of width 5 My (Fig. 9). Each rock uplift history starts at 80 Ma and is constrained not to allow burial (subsidence).

We define a rectilinear grid of nodes that is appropriately oriented with respect to the Alpine fault (Fig. 4), which is the main crustal manifestation of the plate boundary and defines the northwest margin of Fiordland (Barnes et al., 2005). We choose a grid spacing of 50 km in a direction parallel to the Alpine fault, and 25 km perpendicular to the Alpine fault; the choice of distances yields an adequate number of data to fit at every node. For any specific age, the rock uplift value is smoothly interpolated at points other than the grid nodes, specifically datum locations, using natural bicubic splines.

A one-dimensional model is used because the perturbation of isotherms from the horizontal is minor at the relevant depths, as demonstrated by temperature measurements from the Manapouri tunnel (below) and from our previous three-dimensional models incorporating actual Fiordland topography (House et al., 2002). The model thermal history at each datum, T_i(t), is calculated from the exhumation history using a one-dimensional Crank-Nicholson finite difference method with Dirichlet boundary conditions of constant sea-level temperature of 10 °C and a constant tempera ture of 760 °C at 30 km depth. We assume a constant initial geothermal gradient of 25 °C/km, which is consistent with heat-flow measurements from stable locations near Fiordland (Funnell and Allis, 1996); and a thermal diffusivity of 1.0 × 10⁻⁶ m² s⁻¹, which is typical for high-grade metamorphic and igneous rocks similar to those in Fiordland (Clauser and Huenges, 1995; Vosteen and Schellschmidt, 2003).

The upper boundary condition leads to a definition of effective depth, D_i(t), for datum i at age t as the model rock uplift value interpolated at the location of the datum, minus the elevation at the datum. This correction is for the constant temperature experienced during the rock uplift above sea level.

A model thermochronologic age, A_i*, was calculated at each datum from T_i(t). Apatite (U-Th)/He (AHE) and zircon (U-Th)/He (ZHE) ages were calculated using TERRA (Ehlers et al., 2005). Apatite fission track (AFT) and zircon fission track (ZFT) ages were calculated by a code (K. Gallacher, personal commun., 2003) that uses published annealing parameters for apatite (Laslett et al., 1987) and zircon (Tagami et al., 1998). K/Ar data are only relevant (ages within the range of interest) for the most northern node of our model, and there are only a few relevant data available (Claypool et al., 2002). Diffusion kinetic parameters for the white mica K/Ar system are not accurately known; we calculated model ages assuming a constant closure temperature of 350 °C (Dodson, 1973; Hames and Bowring, 1994), though we note that the mylonitic rocks sampled have been subjected to recrystallization at greenschist facies, which could cause us to slightly overestimate the temperature at which the measured age applies (Dunlap, 1997).

For each datum, the interpolated rock uplift history was perturbed by a fixed amount and the thermal history and model ages were recalculated (Fig. 9). A bisection scheme was used to find the value of the rock uplift residual, δr_i, which gives a model age equal to the measured age A_i.

If the cooling age of a datum is older than the maximum age of interest, A_max, then the age relates, at least in part, to events that occurred before the time of interest, and our chosen parameterization of exhumation history is likely to be inappropriate (Fig. 16). However, the datum places a constraint on the cooling history during the period of interest; the sample cannot have exceeded the closure temperature of the respective system during that time.

We treat old data (A_i > A_max) as follows: we redefine each old datum to the maximum age, i.e., A_i = A_max; we determine rock uplift residuals in the usual fashion, giving a model age equal to A_max; the rock uplift residual is only used to discriminate against models that are inconsistent with the datum (see definitions of quality-of-fit metrics below).

In our analysis, we consider the Oligocene–Holocene history (A_max = 34 Ma) as relevant to the process of subduction initiation beneath Fiordland, based upon published analyses of plate boundary motion (Cande and Stock, 2004; Sutherland, 1995) and local studies (Lebrun et al., 2003; Norris and Turnbull, 1993; Turnbull and Uruski, 1993).

Unfortunately, it is not computationally feasible to fully explore the structure of the 36 parameter space that defines the regional uplift model. It would also be a significant challenge to usefully communicate or visualize the results and covariation revealed by any such search. Instead, we consider an alternate simplified approach and explore the geological variability in exhumation history associated with each node independently (below).

We assume a null hypothesis that the optimized model is correct. Our alternate hypothesis is that a new set of 3 parameters that defines the rock uplift history at the node, with the same geological noise, can explain the observations at 5% significance level.

We use a grid search to find, for each node, an envelope of models that are not significantly different from the optimized model.