Tectono‐Thermal Evolution of the Hope‐Kelly Fault System, Southern Alps, New Zealand: Insights From Topographic Analysis and (U‐Th)/He Thermochronology

The fast‐slipping Alpine (∼30 mm/yr), Hope (∼10–20 mm/yr) and Kelly (∼6 mm/yr) faults in the South Island of New Zealand form a complex intersection zone that accommodates tectonic strain along the Australian‐Pacific plate boundary. Analysis of digital topography reveals evidence for stream capture, drainage divide migration, landscape responses to incipient fault development, and preserved enclaves of relic topography that collectively reflect complex interplays between active faulting and landscape evolution. (U‐Th)/He thermochronology of zircon (ZHe) and apatite (AHe) is used to investigate the low‐temperature thermal evolution of rocks in the intersection zone. Weighted mean sample ages for ZHe single grain ages (n = 13 samples) range from ∼9 to 2 Ma, and AHe multi‐grain and single grain aliquot ages (n = 9 samples) range from ∼1.5 to 0.5 Ma. Inverse and forward thermal history modeling reveals distinct spatiotemporal variations in thermal histories. Late Miocene exhumation rates (∼0.6–3.5 km/Myr, assuming geothermal gradients of 33–40 °C/km) through crustal depths of approximately 5–6 km, are interpreted to be controlled by proximity to the Alpine fault, with rocks proximal to the fault recording faster exhumation rates relative to distal samples. Establishment of the Hope‐Kelly fault system in the Quaternary structurally juxtaposed rocks with discordant cooling histories. Rocks throughout the study region record increased cooling rates from ∼2 Ma. Possible causal mechanisms include, spatial changes in rock uplift associated with transport toward the Alpine Fault, increased erosion rates associated with Quaternary climate change, or increased rock mass erodibility associated with development of the Hope‐Kelly fault system.


Orogenic Structure
The Southern Alps formed by deformation and crustal thickening of the Pacific plate along the Pacific-Australian plate boundary ( Figure 1c) and associated uplift and exhumation of Permian to Cretaceous graywacke, cherts and shales of the Aspiring lithologic association and Torlesse composite terrane (Nathan et al., 2002). The Torlesse composite terrane is overthrust and exhumed along the Alpine fault (e.g., Beyssac et al., 2016;Little et al., 2005;Wellman, 1979) and other faults associated with the Southern Alps orogen (Litchfield et al., 2014). Slip on the Alpine fault and associated rock uplift and exhumation rates are generally highest in the central Southern Alps and decrease to the north and south (Beavan et al., 2004;Herman et al., 2010;Little et al., 2005;Michailos et al., 2020). Metamorphosed rocks from the deep crust are exhumed along the Alpine fault, with the highest-grade rocks (amphibolite facies, max temperature ≥640°C; Beyssac et al., 2016;Cox & Barrell, 2007;Nathan et al., 2002) exposed closest to the fault trace. The metamorphic grade decreases systematically eastward to sub-greenschist rocks 15-20 km from the Alpine fault that never reached temperatures greater than 330°C (Beyssac et al., 2016;Cox & Barrell, 2007;Nathan et al., 2002).
Upper-crustal geothermal gradients locally reach up to 125 ± 55°C/km along the Alpine fault (Sutherland et al., 2017) and decrease eastwards to ∼80°C/km near the Main Divide in the highest uplift rate central Southern Alps and ∼30-35°C/km near the Main Divide to the north and south (Michailos et al., 2020;Toy et al., 2010). In contrast to the central Southern Alps (e.g., Shi et al., 1996), isotherms constructed from thermochronologic and seismologic data by Michailos et al. (2020) show that within 20 km south of the Kelly fault there is little to no upward warping of the isotherms with increasing proximity to the Alpine fault ( Figure 1c).
South of the MFS, the Southern Alps are considered to be in large-scale topographic steady state, where rock uplift rates are balanced by erosion rates (Adams, 1980;Batt & Braun, 1997 Cox and Findlay (1995), Little et al. (2005), Michailos et al. (2020), and Warren-Smith et al. (2016). Northern Southern Alps specific aspects include lower (2,000 m) maximum topography and less warping of isotherms compared to geologic cross-sections of the central Southern Alps. Little et al., 2005;Willett & Brandon, 2002). Orographic rainfall and the Alpine fault-driven uplift gradient on the western side of the range maintains a steep, densely forested landscape presently dominated by fluvial incision and landsliding, while the eastern side is characterized by a drier climate, generally lower regional slope, reduced uplift rate, and preservation of glacial U-shaped valleys (Batt & Braun, 1999;Herman et al., 2010;Koons, 1989Koons, , 1990).

Topographic Analysis
We develop and apply two new metrics (Topographic Sinuosity Index, Topographic Integral Index; see Supporting Information S1 for detailed methods) to quantify the topography of a fault trace and its relationship to the surrounding swath topography as a proxy for the evolutionary state of fault-parallel topography. The Topographic Sinuosity Index is the normalized difference between the true length of the fault trace topographic profile and the projected length of the fault trace in map view, describing the number and magnitude of topographic features the fault trace crosses, such as non-parallel ridges or passes. This metric is similar to the commonly used Mountain-Front Sinuosity Index (Keller & Pinter, 1996), although in the vertical rather than horizontal plane. The Topographic Integral Index is the ratio of the area under the curve of the fault topographic profile to the area under the curve of the minimum swath topography and describes how elevated the fault trace is compared to the minimum topography. Low Topographic Sinuosity and low Topographic Integral Indices quantify morphology that may result where enhanced erosion along the fault trace and progressive disruption of drainage networks has produced fault-parallel valleys with the fault trace located near the valley base. To account for uncertainties in Topographic Sinuosity and Topographic Integral measurements, primarily relating to the accuracy and precision of mapped fault traces relative to the topography, we include a ±5% error margin to each measurement.
Because of the systematic difference between the landscape morphology on the western and eastern sides of the orogen, the planform arrangement of watersheds and drainage divides can also provide insights into the status of the landscape response to the Hope-Kelly fault system.

Thermochronology Sampling Strategy
Outcrop samples ( Figure 1b) for (U-Th)/He thermochronology were collected on foot and by helicopter over multiple field seasons between November 2017 and March 2019. The study area is covered with dense vegetation and thick soils that restrict access to unweathered bedrock, so the sampling strategy was partially opportunistic along riverbanks and mountain tops, with little to no rock exposure on the most accessible mountain slopes. Where possible, sample sets were collected at similar elevations and distances from the Alpine fault, and on opposing sides of structures of interest (i.e., Hope and Kelly faults) to examine their effect on rock cooling histories. Rocks were also sampled at the top and base of the highly fractured Sackung Hill and Kelly Range within the Hope-Kelly fault system, to investigate if the dense rock fracturing resulted in detectable increased exhumation rates. Although we sampled in areas affected by faulting, we avoided sampling highly deformed fault rocks or near significant veins to avoid the potential effects of frictional and hydrothermal reheating. Sample lithologies are semi-schists and metagreywackes from the Rakaia terrane of the Torlesse composite terrane (Nathan et al., 2002).

(U-Th)/He Thermochronology
The ZHe and AHe systems are based on the production of 4 He during radioactive alpha decay of 238 U, 235 U, 232 Th, and 147 Sm. However, the decay of 147 Sm generally contributes <1% of total 4 He, compared to typical He age precisions of a few percent (Farley & Stockli, 2002), and is thus considered statistically insignificant when analyzing U and Th-rich minerals, such as zircon. Retention of radiogenic alpha particles ( 4 He) in the crystal lattice is principally a product of time, temperature, crystal size, and accumulated radiation damage (Farley, 2000;Flowers et al., 2009;Guenthner et al., 2013;Reiners et al., 2004). However, a range of additional crystallographic and compositional factors can affect 4 He diffusion and result in overly dispersed ages (Danišik et al., 2017;Gautheron et al., 2009;Gerin et al., 2017;Wildman et al., 2016;Zeitler et al., 2017). Typically, in apatite, 4 He diffusion significantly accelerates at temperatures above ∼40°C and occurs nearly instantaneously at temperatures typically 5 of 27 greater than ∼80°C (Wolf et al., 1998). However, common amounts of accumulated radiation damage in apatite can increase the closure temperature of the AHe system to as much as ∼115° C (Flowers et al., 2022). Within this temperature range, referred to as the He partial retention zone (PRZ), 4 He is partially retained (Spiegel et al., 2009). The higher temperature sensitivity ZHe system has a typical PRZ of ∼130°C-200°C (Biswas et al., 2007;Reiners & Brandon, 2006;Wolfe & Stockli, 2010), but is lower in highly radiation-damaged grains (e.g., Guenthner et al., 2013).
Apatite and zircon grains were separated from rock samples via standard mineral separation and preparation techniques using a jaw crusher, plate mill, and Wilfley table, followed by magnetic and heavy liquids (sodium polytungstate and methylene iodide) separation at the University of Melbourne following the protocol of Kohn et al. (2019). Individual grains were hand-picked and screened for visible inclusions, with euhedral grains preferentially selected where possible. ZHe and AHe data were obtained following the analytical method outlined in Gleadow et al. (2015), except that for ZHe analysis 233 U and 229 Th spikes were used. Further, Fish Canyon Tuff zircons (Gleadow et al., 2015) and Durango apatites (McDowell et al., 2005) were also run as secondary reference materials with each batch of zircons or apatite samples and served as an internal check on analytical accuracy. In this work, the calculation of uncertainties for (U-Th)/He analyses was carried out in two steps. The total analytical uncertainty (TAU) was calculated as the propagated weighted uncertainties on the He, U, Th, and Sm measurements. These were then used to calculate the error on the uncorrected (U-Th)/He ages. The corrected (U-Th)/He age uncertainty was calculated by propagating the TAU value together with an uncertainty of 5% for the uncertainty of the alpha-ejection F T correction value (see Tables 1 and 2).
All ZHe ages reported herein are single grain ages from euhedral zircon crystals with two terminations. ZHe grains were considered outliers if they differed in age from their intra-sample single grain counterparts by more than twice the 2σ age uncertainty. By contrast, low U and Th concentrations, young closure ages, and resulting low radiogenic 4 He content of apatites in this study required analysis of multi-grain aliquots comprising 2-3 apatite crystals in most cases. Most analyzed apatite grains had zero crystal terminations. In the text and figures, ZHe and AHe data are reported as α-ejection corrected ages (after Farley et al., 1996), with intrasample aliquot dispersion described using weighted means with 95% confidence interval (CI) uncertainty (Vermeesch, 2018) and interquartile ranges (IQR), a robust measure of statistical dispersion in He data sets (Boone et al., 2019).

Topographic Analysis Results
Topographic long-profiles of the main traces of the Wairau, Awatere, Clarence, Hope, and Kelly faults and fault-parallel swath topographic minima and maxima (Figures 2a and 2b) reflect the relationships between the fault trace and topography. The Wairau fault is the oldest of the analyzed structures, with the greatest cumulative slip. It has a Topographic Sinuosity of nearly 0.13 ± 0.01 and a Topographic Integral Index of 1.11 ± 0.05 (Figure 2c), indicative of a fault trace that has minimal topographic relief and resides at low elevations in the base of the valley (Figure 2b). The Awatere and Clarence faults have intermediate values of both metrics (Figure 2c); the deviation from "mature" values is primarily due to interaction with the Inland and Seaward Kaikoura Ranges on the east side of the island. The topography in these areas is still primarily fault-parallel and fault-controlled (Collett et al., 2019;Duvall et al., 2020). The Hope fault has a Topographic Integral Index value (1.17 ± 0.05), slightly higher than the Awatere fault (1.20 ± 0.06), and a Topographic Sinuosity Index value (0.57 ± 0.03) close to that of the Wairau fault (0.13 ± 0.01; Figure 2c). This indicates that the fault trace crosses relatively few topographic undulations but is located at a higher elevation than local minimum topography. Visual inspection confirms this is primarily due to interaction with the Kaikoura Ranges at the Hope fault's eastern end (Figure 2b).
In contrast to the other faults, the Kelly fault has the highest Topographic Sinuosity Index (3.82 ± 0.19) and the highest Topographic Integral Index (1.41 ± 0.07), suggesting a topographically immature fault (Figure 2c). The Kelly fault trace transects some of the highest topography within the swath area (high Topographic Integral value) and traverses multiple fault-perpendicular ridges and valleys (high Topographic Sinuosity; Figure 2b). Ridges within the Hope-Kelly fault system exhibit dense networks of faults that reflect a combination of gravitational collapse and tectonism (Beck, 1968;Pere, 2009;Vermeer et al., 2021). The absence of Kelly fault-parallel topography and the prevalence of dense fault networks suggests that the Kelly fault is incipient and the topographic response to faulting has been insufficient to erode pre-existing topography. The landscape surrounding the Kelly Weighted mean age ± 95% confidence interval (Ma) (MSWD) 28.9 ± 1.1 (1.45) Note.
Italicized ages were not included in weighted mean age calculations, which were calculated using IsoplotR (Vermeesch, 2018 Gleadow et al., 2015). Where more than one grain has been analyzed then the mass-weighted average radius (MWAR) of apatite crystals measured in the aliquot analyzed is given. g Standard deviation of the MWAR is used as a guide for the "tightness" of the range of single crystal radii picked from more than one grain within a sample. h Grain morphology -0T = no terminations, 1T = one termination, 2T = 2 terminations. i AHe age corrected for alpha ejection after Farley et al. (1996) Planform analysis of topography around the Hope-Kelly fault system is also indicative of a disrupted and dynamic landscape (Figure 2d). In the central Southern Alps, where a regional steady state has been achieved, the main divide position is stable relative to the Alpine fault (He et al., 2021;Herman & Braun, 2006;Willett & Brandon, 2002). The central Southern Alps has distinct drainage patterns on either side of the Main Divide; the west side has primarily fluvial dendritic E-W trending drainage systems, while the east side has glacial, N-S trending valleys, with the main divide separating the two domains (Craw et al., 2013). The steeper rivers and higher erosion rates on the west side counteract the westward advection of rocks, and as a result, the main divide is primarily erosional and at the highest elevation (Figure 2e; He et al., 2021). East side features may be occasionally captured and incorporated into the west side drainages (i.e., Perth River or Landsborough River, Craw et al., 2003Craw et al., , 2013Herman & Braun, 2006), forming a wind gap or relatively low elevation pass across the divide (Figure 2e). However, these capture events are likely isolated to single drainages within the central Southern Alps, as analog and digital modeling of similar systems have shown that the interplay of erosion, uplift, and advection will maintain large-scale topographic balance and main divide position relative to the Alpine fault surface trace until perturbed (He et al., 2021;Herman & Braun, 2006). Within the western MFS, the drainage pattern has a complex overprinting of MFS fault-controlled ENE trending valleys capturing and dissecting pre-existing N-S trending eastern Southern Alps style valleys (Craw et al., 2013). Just south of the Hope-Kelly fault system, the Main Divide has a ∼20 km eastward jog (Figure 2d) that includes many passes (Figure 2e) where NS trending valleys are cut by the main divide, resulting in wind gaps and expansion of the western catchment area. The ∼90° bends in the Taipo and Arahura rivers and the sharp junction of the Otira and Taramakau rivers may be relict main divide topography, where the old EW/NS valley pattern is preserved but captured by west-flowing rivers (Figure 2d). Compared to the steady-state landscape to the south, the relict divide topography is closer to the Alpine fault surface trace (Figure 2d), likely reflecting the lack of headward erosion to balance westward advection of topography with the rocks. These landscape anomalies suggest that this area is currently undergoing a period of disequilibrium and adjustment brought on by the formation of the Hope-Kelly fault system, still in the early stage of imparting its influence on the topography of this region via landscape responses to faulting.

(U-Th)/He Thermochronology Results
Detailed ZHe and AHe results are reported in Tables 1 and 2. Sample weighted mean ages and location information are in Table 3, and are presented geospatially in Figures 3a and 3b, respectively. Detailed sample and (U-Th)/He data are publicly accessible via AusGeochem , where the data can be interrogated geospatially; see https://app.ausgeochem.org/doi/10.58024/AGUM98746938#/.
Single grain ZHe ages (n = 72) from 14 samples in the Hope-Kelly fault region range from 1.00 ± 0.1 to 33.8 ± 1.7 Ma, with weighted mean ZHe ages of between 1.16 ± 1.0 and 7.91 ± 0.5 Ma. In general, intrasample ZHe age dispersion shows no discernible relationship to sample elevation (Figure 4a) or grain size (Supporting Information S1). Single grain ZHe ages do, however, show a negative correlation with effective uranium content (eU; eU = U ppm + 0.235 · Th ppm), a proxy for accumulated radiation damage (Flowers et al., 2022), which ranges from 126 to 3,569 ppm (Figure 4), indicative of increased 4 He diffusion due to accumulated radiation damage (Guenthner et al., 2013). Regionally, ZHe ages increase with distance from the Alpine fault ( Figure 5).
AHe analyses include two single grain ages (JT4: 1.18 ± 0.07 Ma; J115: 1.05 ± 0.07 Ma), and 45 multigrain aliquot ages ranging between 55.1 ± 3.3 and 0.47 ± 0.03 Ma. However, a total of 40 single or multi-grain AHe aliquots yield ages younger than 2.5 Ma, while only three are between 2.5 and 10 Ma and another four aliquots are greater than 10 Ma. For aliquots younger than 2.5 Ma, which comprise the bulk of AHe data, the intrasample age dispersion is low to moderate (IQR values from 0 to 0.75 Ma). Potential sources of uncertainty in AHe ages include grain morphology and He ejection correction (most grains had no terminations and grain morphology is averaged within each aliquot), compositional zoning of parent isotopes, microscopic U-and/or Th-bearing micro-inclusions; external He implantation, as well as defects, vacancies and amount of radiation damage in the crystal structure. Outlier aliquots were designated as those which return an age more than 2x the age of the other intra-sample aliquots or that are older than the ZHe ages from the same rock; under these criteria, all aliquots with AHe ages >3 Ma are considered outliers. There is no apparent relationship between AHe age and eU concentration (6-413 ppm, Figure 4d) or grain size (Supporting Information S1). AHe sample mean ages exhibit a slight  positive correlation with elevation ( Figure 4b) as would be expected for denudational exhumation, although the increase in ages fall within mean age uncertainties. In contrast to ZHe data, AHe multi-grain aliquot and single grain ages exhibit no apparent relationship to distance from the Alpine fault ( Figure 5).

Thermal History Modeling
To quantify the upper crustal cooling history of the study area, inverse (Figures 6 and 7, Supporting Information S1) and forward (Supporting Information S1) thermal history modeling of low-temperature thermochronology data was performed using QTQt version 5.8 (Gallagher, 2012) following the modeling protocol described in the Supporting Information S1. Only samples JT4 and J115 yielded single grain AHe data which could be modeled jointly with ZHe results. While multi-grain aliquot AHe ages were unable to be modeled, time-temperature windows corresponding to intra-sample aliquot IQRs and a conservative ∼40°C-80°C temperature range were plotted onto inverse thermal history modeling results to allow for comparison of time-temperature reconstructions and observed AHe data.

Individual Sample Inverse Thermal History Models
Inverse thermal history modeling was performed via the protocol described in Supporting Information S1. In general, inverse thermal history models reproduce observed (U-Th)/He ages within uncertainty limits, with the notable exception of single grain ages from J135 which exhibits significantly more intra-grain age dispersion (4.1-24.9 Ma), corresponding to a wide range in eU content (2210-207 ppm). This reflects the now well-reported inability of the applied zircon radiation damage accumulation and annealing models of Ginster et al. (2019) and Guenthner et al. (2013) to reproduce highly disperse ZHe data sets (e.g., Boone et al., 2018; Guenthner  Figure 6). Prior to late Miocene cooling, inverse models exhibit widely distributed time-temperature envelopes spanning between 300°C and near-surface temperatures that predate the oldest observed ZHe ages, which are poorly constrained by the data.
The timing of late Miocene cooling through the ZHe PRZ varies systematically in relation to sample distance from the Alpine fault. Thermal history models suggest samples further to the east cooled below ∼220°C as early as ∼10-8 Ma, while samples closer to the Alpine fault cooled through the upper ZHe PRZ as recently as 2 Ma. Inverse thermal history models generated using only ZHe data exhibit monotonic late Miocene-recent cooling paths. However, the Bayesian approach of the QTQt thermal history modeling software favors simpler, prolonged cooling paths (Gallagher, 2012), allowing for the possibility that ZHe data are consistent with more complex late Miocene-recent cooling histories.
By contrast, samples which yielded single grain AHe data that could be modeled (samples JT4 and J115) produced time-temperature reconstructions characterized by a more complex, multi-stage cooling history ( Figure 6). These models instead exhibit an initial late Miocene phase through the ZHe PRZ, followed by a period of thermal stability or slight reheating, before undergoing a renewed cooling phase through the AHe PRZ to near-surface temperatures beginning around 2 Ma.  Herman et al. (2007). Weighted mean ages and 95% CI are shown for each sample. A trendline fit through the weighted means shows a regional trend of older ZHe ages with distance eastwards from the Alpine fault. For AHe the mean of the weighted mean ages passes through the 95% weighted mean confidence interval of most samples. The white line is the maximum mode path. Yellow boxes illustrate the range in AHe dates, where available, corresponding to the AHe IQR and a temperature range of 40°C-80°C. The small graph accompanying each inverse model shows a plot of the observed age versus predicted age for each modeled grain. Inverse models are presented individually in Supporting Information S1.
In all cases, AHe multi-grain aliquot age ranges overlap inversely predicted time-temperature envelopes ( Figure 6).
Despite differences in cooling paths (multi-phased vs. monotonic) recorded in inverse modeling results of samples with and without single grain AHe data, forward thermal history modeling (Supporting Information S1) shows that ZHe data from nearly all samples across the study area are consistent with a more complex multi-staged thermal history involving a cooling hiatus at ∼3 Ma followed by a transition to rapid cooling ∼ 2 Ma (Supporting Information S1).

Grouped Inverse Thermal History Models: Sackung Hill, Kelly Range Tops, and South of Kelly Fault (SKF)
Grouped inverse thermal history models for samples atop the Kelly Range (J98, J105, 1366), samples on the south side of the Kelly fault (SKF: J127 and JT4), and samples from Sackung Hill (J115, J116, J30, J12, J135) exhibit three distinct cooling histories for these regions (Figure 7) that are consistent with their constituent individual time-temperature reconstructions. Samples from Sackung Hill and SKF both exhibit a multi-stage cooling history involving relatively prolonged late Miocene cooling through the ZHe PRZ (∼200°C-130°C; Flowers et al., 2022), followed by a period of slowed cooling between ∼3 and 2-1.5 Ma and a subsequent phase of renewed cooling through the AHe PRZ (∼80°C-40°C; Flowers et al., 2022;Wolf et al., 1998) to surface temperatures. The Kelly Range samples, by contrast, record significantly greater Quaternary cooling than the other two groups, cooling from temperatures above 200°C within the last 2-3 Ma.
Using the expected and lower 95% confidence interval time-temperature paths of grouped inverse models, exhumation rates were calculated based on the modeled isotherms of Michailos et al. (2020), which predict an increase in geothermal gradient of 2-8°C/km in the uppermost 2-3 km of the crust, and assuming a constant paleo-geothermal regime over the last ∼6 Ma (Table 4). . Like the individual inverse models, grouped thermal history simulations are able to reproduce all but the oldest single grain ZHe ages, particularly from sample J135, reflecting the applied zircon radiation damage accumulation and annealing models' inability to reproduce highly disperse ZHe data sets. This model indicates cooling through the HePRZ ∼8 Ma, and an increased cooling rate after ∼2 Ma. The timing of cooling through 200°C and 100°C are shown in thin dashed lines, which are used for estimating exhumation rates.

Late Miocene-Pliocene Exhumation of the Northern Southern Alps
ZHe thermochronology from the northern Southern Alps indicates late Miocene-Pliocene cooling that we interpret to record orogenic rock uplift and associated exhumation. Thermal history modeling indicates that rocks now outcropping in the eastern Hope-Kelly fault region experienced cooling through the ZHe PRZ (∼200°C-130°C) as early as 5-8 Ma (Figures 6 and 7, Sackung Hill samples), indicating an exhumation rate of 0.6-1.1 km/Ma (Table 4, SH sample set). Further west, ZHe data record cooling from similar mid-crustal temperatures over just the last 2 Ma, suggesting exhumation rates >2.5 km/Ma (Table 4 KR sample set), that we attribute to increased rates of exhumation associated with more rapid inferred rock uplift rates above the Alpine fault ramp, and locally increased denudation rates.
A similar east-to-west increase in exhumational cooling rates has been observed in the central Southern Alps as recorded by the zircon fission track system (Herman et al., 2009;Tippett & Kamp, 1993), which has a slightly higher temperature sensitivity range (closure temperature in the range of ∼200°C-260°C; Reiners & Brandon. 2006) than that of the ZHe system. The ZHe mean ages here are slightly older at the eastern extent (up to ∼8 Ma weighted mean age, sample J116, ∼15 km from the Alpine fault) than zircon fission track ages at equivalent distances from the Alpine fault further to the south, a pattern we interpret as reflecting a general northwards decrease in exhumation rates from maxima at Aoraki/Mt Cook, as recognized by prior workers (Little et al., 2005;Michailos et al., 2020;Tippett & Kamp, 1993). While rapid exhumation can advect heat vertically, increasing the geothermal gradient (Koons, 1987;Shi et al., 1996;Warren-Smith et al., 2016), in the northern Southern Alps 3D isothermal surfaces inferred from thermochronology and seismology (Michailos et al., 2020) suggest that little upward warping of low-temperature (<300°C) isotherms has occurred within 5-20 km of the Alpine fault. Based on the slower late Miocene cooling rate, presently exposed rocks located further east of the Alpine fault appear to have exhumed through ZHe closure temperatures whilst located in the eastern sector of the orogen, in a tectono-thermal regime dominated by lateral tectonic transport of rocks, lower uplift rates, and slower regional denudation.

Quaternary Tectono-Thermal Evolution
Despite clear differences in late Miocene cooling histories across the northern Southern Alps related to westward advection of rocks across and up the Alpine fault ramp, upper crustal thermal histories across the study area appear to be consistent in the Quaternary. This is reflected in the lack of a clear east-to-west trend in AHe ages across the western side of the orogen, in contrast to distinct spatial trends in their higher temperature-sensitive Note. Geothermal gradients are extracted from the isotherm of Michailos et al. (2020). The time and temperature ranges are taken from the grouped inverse time-temperature thermal history models in Figure 7. The exhumation rate is calculated by the temperature change divided by the time over which that change took place, divided by the appropriate geothermal gradient. Each exhumation rate range takes into account the range of time and geothermal gradient for each temperature range.
ZHe counterparts ( Figure 5). Joint inverse modeling of ZHe and AHe data (samples JT4 and J115 and grouped models, Figures 6 and 7) and forward thermal history simulations (Supporting Information S1) are consistent with spatially invariant time-temperature histories from approximately 2 Ma. There are a few possible mechanisms by which these observations may be explained.
Structurally controlled spatial variations in the vertical component of rock trajectories through orogenic isotherms could impart significant changes in rock cooling rates as they are laterally advected through the orogen (and by proxy, ZHe and AHe thermochronologic ages; Batt & Brandon, 2002). Under simplified conditions where westward-increasing tectonic rock uplift is balanced by erosion, rocks presently exposed further east (where rock uplift through the ZHe PRZ was driven by relatively slower slip rate faults; Litchfield et al., 2014) could be expected to yield older ZHe thermochronologic ages relative to presently exposed rocks that were exhumed through the ZHe PRZ by the faster slipping Alpine Fault. The potential for steeper geothermal gradients associated with more rapid rock uplift in the Alpine Fault region would only be expected to further increase this ZHe age variance (in a potential trade-off with ZHe PRZ depth and associated thermal gradient). However, progressive westward rock transport of eastern orogenic rocks toward the Alpine Fault could bring rocks that previously cooled slowly through the ZHe PRZ in the eastern part of the orogen into the Alpine Fault dominated tectonothermal regime prior to (or during) cooling through AHe PRZ. Thus, rocks with discrepant ZHe ages (associated with their distinct paleo-tectonic positions during ZHe PRZ-related cooling) could have similar AHe ages related to similar tectonic positions (e.g., within the Alpine Fault-dominated regime) by the Quaternary. In this hypothesis, rock transport westward through the orogen could potentially create the observed ZHe and AHe age distributions without requiring a temporal change of orogen dynamics. However, the AHe ages would still be expected to exhibit an eastward-increasing age distribution (reflecting the westward-increasing gradient of uplift and erosion rates, which peak near the Alpine fault), which is not observed in our data. Finally, thermochronometers with lower closure temperatures are less affected by these orogen dynamics and are more affected by local conditions at the surface (Batt & Brandon, 2002), and understanding exactly how well this mechanism explains the observed ZHe and AHe ages would require structural-thermal modeling beyond the scope of this study.
Another possible hypothesis to explain the spatially invariant AHe ages is that increases in exhumation rates in the Quaternary extended the spatial extent of rapid cooling beyond the near-Alpine Fault region, such that any east-west (tectonic) gradients in rock exhumation rates were reduced to below the resolution of the AHe data. Some previous thermochronological studies propose a widespread increase in Southern Alps exhumation rates around 2 Ma, attributed to the onset of Quaternary glacial cycles and related increased erosion rates (Jiao et al., 2017;Tippett & Kamp, 1993). Others find no evidence for accelerated Quaternary erosion or increased exhumation rate related to glaciation of the Southern Alps (Herman et al., 2007(Herman et al., , 2010Lang et al., 2018Lang et al., , 2020. Increased Quaternary climatic instability (e.g., more rapid glacial-interglacial fluctuations) and associated landscape instability (Peizhen et al., 2001) could have contributed to the apparent increases in rock exhumation rate we observe in this study. Climatically enhanced erosion rates could have affected exhumation rates and AHe ages either independently of, or in concert with, the progressive tectonothermal changes described above, and potential rock strength changes discussed below.
A final possible hypothesis for exhumation rate changes in the study area is a local increase in rock erodibility (i.e., decrease in rock mass strength) associated with the development of the Hope-Kelly fault system. The Hope fault is estimated to have initiated between 1.6 and 0.6 Ma based on Holocene slip rates applied to displaced bedrock units east of the Main Divide (Langridge et al., 2013). As the Hope-Kelly fault system formed and accumulated displacement, built-up rock fracturing and increased frequency of strong shaking may have enhanced local erosion rates both in the principal slip zones of the faults and in the fractured surrounding rock mass, thereby locally increasing denudation rates, exhumation, and rock cooling rates. The progressive development of the structurally complex interaction zone between the Hope, Kelly and Alpine Faults could have shifted the region of faulting-affected rock mass weakening from the near-Alpine Fault region (prior to Hope-Kelly Fault establishment) eastward (to the region of the present interaction zone). This could have allowed surface processes to destabilize rock masses more effectively, via changes to rock mass strength and rearrangements of fluvial geomorphic networks.
Distinguishing which of these hypotheses played the most prominent role(s) in generating the observed structural-thermochronological characteristics identified herein is beyond the resolution of our current data sets. Further modeling of our existing results is limited by the number of presently acquired single grain AHe ages.
However, the application of other complimentary low-temperature thermochronometers, such as apatite fission track analysis, or numerical geodynamic-landscape evolution modeling techniques may provide further insights into the geological mechanisms governing the upper crustal tectono-thermal evolution of the greater Hope-Kelly fault region.

Total Displacement of the Hope-Kelly Fault System
The total displacement of the Hope fault is 8-13 km (Langridge & Berryman, 2005;Nathan et al., 2002), based on bedrock displacements on the central Hope fault, east of the Main Divide and Hope-Kelly-Alpine intersection zone. Apparent dextral displacements of the metamorphic zones (Figure 1b) within the Hope-Kelly fault system are more difficult to relate directly to the age and total displacement of the Hope or Kelly faults because of the oblique fault kinematics and interaction with the Alpine fault ramp. The >5 km apparent dextral offset of the east-dipping metamorphic zones across the westernmost Hope fault was produced by normal kinematics on the E-W striking fault, while the <2 km of apparent offset of the metamorphic zones across north dipping oblique normal-dextral Kelly fault splays may underestimate net total displacement (Nathan et al., 2002;Vermeer et al., 2021). The juxtaposition of young (1-2 Ma) and older (3-4 Ma) ZHe ages from the Kelly Range and adjacent region on the south side of the Kelly fault (SKF group) also provide first-order constraints on the total displacement on the central part of the Kelly fault. Retro-deforming the slip on the Kelly fault restores the SKF samples to positions further to the east and more distal from the Alpine fault relative to the Kelly Range samples (Figures 3 and 8). The contrasting paleo-positions of these presently proximal sample groups can account for these crustal blocks having experienced distinct cooling histories (Figure 8a). Applying the late Pleistocene-Holocene slip-rate of the Kelly fault of ∼6 mm/yr (Vermeer et al., 2022) to back slip on the Kelly fault for 2 Myr yields a maximum of 12 km of dextral displacement, sufficient separation to have produced the distinct cooling histories of the Kelly Range and SKF groups. The slower late Miocene (ZHe cooling through the PRZ ∼8 Ma) cooling recorded by Sackung Hill samples may provide broad longitudinal paleo-position constraints for SKF samples JT4 and J127 prior to the formation of the Hope-Kelly fault system, which record faster cooling (ZHe cooling through the PRZ between ∼8 and 6 Ma) and suggest that prior to Kelly fault formation these samples were on a trajectory to surface further west than Sackung Hill samples and further east than the Kelly range samples. This suggests that the total Kelly fault dextral displacement is likely less than 20 km.

Topographic Response to the Hope-Kelly Fault System
Topographic analysis indicates that the landscape is dynamically evolving in response to Hope-Kelly fault system development. The southern Kelly fault splay has elevated topographic indices compared with the other MFS faults (Figure 2c), including the westernmost Hope fault. The westernmost section of the Hope fault, where the strike is ∼90°, has a well-developed fault-parallel valley (Taramakau River), but unlike the rest of the Hope fault system and other major MFS faults, it is primarily normal slip (Vermeer et al., 2021). This exemplifies how fault kinematics may promote evolution and maintenance of fault associated large-scale topography. This also shows that care must be taken in applying these fault trace topographic metrics to faults with various kinematics, as the resulting indices likely have different significance for strike-slip, normal and reverse faults.
In the eastern part of the Hope-Kelly fault system, the Taramakau River valley trends ∼070 and is sub-parallel to the Hope and Kelly faults, an orientation that is not relict and suggests its formation has been encouraged by the presence of the faults. However, the Kelly fault west of the Otehake River (Figures 2b and 2d) has not yet succeeded in forming a fault-parallel valley. Prior to the Hope-Kelly fault system's disruption of the erosion-advection balance, the Main Divide and topographic maxima were probably located ∼15-20 km from the Alpine fault surface trace, where Sackung hill is currently (Figure 2d). Dense fracture networks on Sackung Hill, the Kelly Range and nearby hills are likely a sign that this topographic arrangement is not stable with the relatively young fault system and they are in the process of being broken down and may eventually be replaced by a fault-parallel valley. This Hope-Kelly fault driven increase in erosion rates, especially near the Main Divide (formation of the along-fault upper Taramakau River valley and capture of east-flowing rivers and preceding erosion of the drainage divide), may be the topographic signature of processes which could decrease the previous westward-increase of erosion rates and resulted in the observed spatial consistency of AHe ages (Section 6.2).
The topographic response to increased erosion and discrete displacements of the Hope-Kelly fault system will also compete with the uplift, advection, and erosion controlled by the Alpine fault itself. This may be especially important west of the Otira River where the Kelly fault begins to splay (Vermeer et al., 2021), reducing the incremental slip that is localized on any one fault trace, thus reducing localized erosion and drainage deflection, and distributing fracturing and erodibility across the landscape. Meanwhile, Kelly fault parallel erosion processes must compete with westward increasing erosion and uplift rates driven by the Alpine fault and the orographic rainfall. These factors may slow the development of Kelly fault-parallel topography, to the extent that the distributed fault zone may never sufficiently concentrate erodibility onto one or more fault splays, and a continuous Kelly fault parallel topographic signature might never form within the western Southern Alps. This implies that regardless of age, the western Kelly fault(s) may never form a characteristic fault parallel valley, and if so, it would never be recognized by mapping techniques limited to large-scale topography. Thus, mapping of bedrock deformation and localized fault scarps may be required to recognize a fundamentally distributed fault system, even as it matures and accumulates tens of kilometers of slip.
Although the southern Kelly fault splay does not have a large-scale topographic signature, the Hope-Kelly fault system has managed to disrupt the advection-erosion balance of the Main Divide, promoting capture of portions of the N-S trending Waimakariri River headwaters by the west flowing Taramakau, Taipo, and Arahura rivers (Figure 2d). The combination of ∼N-S extension of the Hope-Kelly fault system and headward erosion of the Taramakau River along the narrow eastern end of the Hope-Kelly fault system are potential mechanisms for Main Divide eastward migration and river capture. Although the Kelly fault splays have not yet, and may never develop individual topographic signatures, the fault system overall has disrupted the topographic balance of the orogen.

Conclusions
Late Miocene to Pliocene exhumation recorded by ZHe thermochronology exhibits cooling rate variations relative to the Alpine fault, reflecting its fundamental role in the exhumation history of the northern Southern Alps prior to the Quaternary. In contrast, AHe data exhibit no apparent age trends with respect to distance from the Alpine fault, position relative to the Hope-Kelly fault system, or elevation. AHe sample weighted mean ages and thermal history modeling indicate increased cooling rates in the Quaternary, but exactly how this relates to concurrent tectonic, structural and climatic processes which may have affected denudation rates during orogenesis remains ambiguous. Lateral variation in the vertical component of rock trajectories associated with tectonic transport from the eastern part of the Southern Alps orogen into the western Alpine fault ramp-associated zone of enhanced rock uplift provides a feasible mechanism to explain the observations made in this study. However, we cannot refute the possibility that climatically driven landscape instability and reductions in rock mass strength associated with establishment of the Hope-Kelly fault system could have further contributed to apparent increases in Quaternary rock exhumation rates and overprinting of the expected spatial trend in AHe ages. Re-arrangements in the geometry of fluvial systems in response to fault system developments may have further contributed to rapid exhumation rates. Distinguishing between the complex interplay of tectonic, climatic, and rock strength contributions to the evolution of orogenic landscapes remains an important challenge here and in other orogens worldwide.

Data Availability Statement
All data (AHe and ZHe analyses) used in this study are available in the data tables and at AusGeochem, https:// app.ausgeochem.org/doi/10.58024/AGUM98746938#/ under the Data Package Vermeer et al. (2023) (https://doi. org/10.58024/agum98746938). All software used in the analysis is available from the original sources, as cited in the text and references.