The Significance of Heat Transport by Shallow Fluid Flow at an Active Plate Boundary: The Southern Alps, New Zealand

Fluid flow can influence fault behavior. Here we quantify the role of groundwater heat advection in establishing the thermal structure of the Alpine Fault, a major tectonic boundary in southern New Zealand that accommodates most of the motion between the Australian and Pacific Plates. Convergence on the Alpine Fault has rapidly uplifted the Southern Alps, resulting in high geothermal gradients and a thin seismogenic zone. A new equilibrium temperature profile from the 818‐m‐deep Deep Fault Drilling Project 2B borehole has been interrogated using one‐dimensional analytical models of fluid and rock advection. Models indicate a total heat flux of 720‐mW m−2 results from groundwater flow with Darcy velocities approximating to 7.8 × 10−10 m s−1. Groundwaters advect significantly more heat than rock advection in the shallow orogen (<6‐km depth) and are the major control on the subsurface temperature field.


Introduction
Fluid flow can influence the nature and timing of failure on faults by changing temperatures and pore-fluid pressures, and by modifying the mechanical properties of fault materials through chemical alteration and mineralization [Hubbert and Rubey, 1959;Sibson, 1973Sibson, , 1990Kohlstedt et al., 1995]. The Alpine Fault is a major plate-bounding fault accommodating ~75% of the ~10 mm·yr -1 convergent and 37 mm·yr -1 dextral strike-slip motion between the Australian and Pacific plates in South Island, New Zealand (Fig. 1a) Cooper, 2001, 2007;DeMets et al., 2010]. Although the relative plate motions are high, no large (M w > 7) earthquakes have occurred on the central Alpine Fault since European settlement ca. 1800 AD, and there is no evidence for aseismic creep at the surface [Sutherland et al., 2007]. However, paleoseismic records indicate that the Alpine Fault has ruptured in multiple large (M w > 7) to great (M w ≥ 8) earthquakes every 291 ± 23 years with the most recent event in 1717 AD [Wells et al., 1999;Sutherland et al., 2007;Berryman et al., 2012;Howarth et al., 2016;Cochran et al., 2017]. Consequently, the Alpine Fault appears to be late in its seismic cycle. Convergence across the Alpine Fault has resulted in rapid uplift of the >3500 m Southern Alps and earlier models, supported by micro-seismicity observations, have predicted high geothermal gradients and the uplift of the brittle-ductile transition to shallow levels, and consequent thinning of the seismogenic crust [Allis et al., 1979;Koons, 1987;Leitner et al., 2001;Townend et al., 2012]. Here we interrogate a new equilibrium temperature profile from the 818 m-deep International Continental Scientific Drilling Program (ICDP) Deep Fault Drilling Project 2B (DFDP-2B) borehole to quantify the relative roles of heat conduction, heat transfer by rock advection (defined here as movement of rock and entrained pore water relative to a fixed reference frame), and heat transfer by fluid flow (movement of water relative to rock) in establishing the local, near-surface geothermal profile of a major orogen-bounding, continental strike-slip fault.
The ~45 degree-dipping Alpine Fault has uplifted and juxtaposed amphibolite facies quartzofeldspathic rocks of the Alpine Schists against Cretaceous granitoids and Paleozoic metasediments in its footwall [Cox and Barrell, 2007], with vertical uplift rates of 6-9 mm•yr -1 [Little et al., 2005]. Heavy orographic precipitation (~10 m•yr -1 ) on the western side of the Southern Alps results in high rates of physical erosion and has produced a quasi-steady-state topographic regime in which the range of the highest topography, the Main Divide, runs parallel to but offset 15-20 km southeast from the Alpine Fault, connected by a series of valleys and 1000 to 1500 m-high ridges running normal to the fault (Fig. 1b) [Koons, 1989[Koons, , 1990.
A >50 m thick zone of intense hydrothermal alteration in the hanging wall of the Alpine Fault zone [Cooper and Norris, 1994;Warr and Cox, 2001;Sutherland et al., 2012] and numerous warm (<60°C) geothermal springs in the Southern Alps ( Fig. 1b) [Barnes et al., 1978;Menzies et al., 2014;Cox et al., 2015] issuing meteoric water attest to significant fluid flow penetrating to depths that may influence orogenic processes. The topography of the Southern Alps drives two significant components of groundwater flow (Fig. 1c): (1) flow from ridges into valleys that produces a relatively large upward vertical component to flow beneath the valleys; and (2) flow parallel to ridges and valleys, down the topographic gradient from the Main Divide at ~3500 m to the Alpine Fault trace at ~100-300 m elevation. Differences in hydrological heads across the fault slip zones in the DFDP-1 boreholes [Sutherland et al., 2012], isotopic measurements [Menzies et al., 2016], and low permeabilities measured for the principal slip zone materials [Boulton et al., 2012;Carpenter et al., 2014] indicate that the Alpine Fault is a barrier to flow and likely separates the Southern Alps groundwater from that of the footwall plain. Consequently, the down-valley component is driven upward upon reaching the fault near valley mouths.
The DFDP-2B borehole is located near the mouth of the Whataroa Valley (Fig. 1b). The rock sequence sampled in the borehole consists of (from base to top): Alpine Schistderived mylonites, protomylonites and schist (Layer A); lacustrine diamictite and silty sediments (Layer B) and alluvial gravels (Layer C, Fig. 1d, 2). The borehole has a curvature toward the fault. An equilibrium temperature profile of the borehole, with measurements at 1 m intervals, was obtained using an optical fiber down its length, following repeated measurements during relaxation of the drilling induced perturbation [Sutherland et al., 2017] (see Supporting Information for further information). As the borehole was closed to flow by cemented casing along its entire length during collection of the temperature data, all of which were collected >6 months after drilling, we assume that observations are representative of the background state and are unaffected by the borehole. In the absence of creep strain on the Alpine Fault [Sutherland et al., 2007], we consider observation points as fixed in space. The profile has an average geothermal gradient of 125±55 °C•km -1 [Sutherland et al., 2017], with three distinct discontinuities in gradient separating the profile into four sections (Fig. 2a). The shallower two of these discontinuities correlate with changes from Alpine Schist rocks to valley-filling sediments and from lacustrine sediments to alluvial gravels. The third discontinuity corresponds to a spike in measured radon gas during drilling [Townend et al., 2017] but no discernible change in rock type. Conditions proximal to the borehole are of significant vertical fluid flow due to both components of the regional groundwater flow regime, convergence into the valley and against the fault. This was confirmed by rapid increases in borehole hydraulic head with depth, rising to >60 m head at 818 m [Sutherland et al., 2017].  Cooper, 2001, 2007;DeMets et al., 2010]. b. The Whataroa Valley is one of a series of valleys running perpendicular to the Alpine Fault in the central Southern Alps, bounded by the highest topography at the Main Divide. Many contain warm springs [Barnes et al., 1978;Menzies et al., 2014;Cox et al., 2015]. The DFDP-2 boreholes are located at the Whataroa Valley mouth, with boreholes drilled during the first phase of the project, DFDP-1, in nearby Gaunt Creek. c. Topographic flow from ridges into valleys and away from the Main Divide, driven up along the impermeable Alpine Fault contribute to up-flow at valley mouths. d. Illustration of the rock sequence penetrated by the DFDP-2B borehole (see Toy et al., [2017] for further detail). Projection of units away from the borehole is intended only to illustrate the geological setting.

1D Thermal Modelling and Analysis
Past thermal models for the Southern Alps considered only heat transport through conduction and rock advection [Allis et al., 1979;Koons, 1987;Allis and Shi, 1995]. Recent regional 3D models [Sutherland et al., 2017] also incorporate heat transport by groundwater circulation and are able to simulate high geothermal gradients beneath the Whataroa Valley, similar to those observed at DFDP-2B. We build upon this earlier work to obtain estimates for fluid and heat fluxes at depth, to quantify the importance of fluid flow as a means of heat transfer in the shallow Alpine Fault zone and to provide insights into the local fault zone structure and fluid flow regime. We fit the DFDP-2B temperature profile using the simplest suitable models, stripping away complexity to gain insights into the hydrothermal processes influencing the Alpine Fault. On the basis of a priori knowledge of the system, indicating steep vertical temperature gradients in the Southern Alps and regional converging fluid flow beneath DFDP-2B, we assume that the conductive heat flux is close to vertical for the measured profile. Implicit in this assumption is that temperature variations associated with borehole curvature are negligible. Variations in thermal diffusivity have the potential to produce discontinuities in temperature gradients. However, the most prominent change in DFDP-2B thermal gradient at 691 m occurs without any change in rock type and would require an unrealistically large shift in thermal diffusivity to explain by rock thermal properties alone (see Supporting Information). We therefore investigate the potential of fluid flux changes to produce these discontinuities, in the absence of changes in rock properties. In this approach, discrete kink points within the data represent locations where lateral flow pathways significantly modify the fluid and heat flow up the profile. We do not rule out some influence of changes in thermal properties on the data and we consider possible properties for rocks sampled in DFDP-2B in our uncertainty range for fixed rock properties [Vosteen and Schellschmidt, 2003;Mottaghy et al., 2005Mottaghy et al., , 2008.
On this basis, we present one-dimensional models for fluid and heat flow. Our approach is similar to that used by Bredehoeft and Papadopulos [1965] for estimating fluid flow rates, and that used by Allis et al. [1979] for modelling the thermal consequences of rock uplift, but incorporates both processes. Although the approach of Bredehoeft and Papadopulos [1965] can be modified to include a term representing lateral heat transfer [Lu and Ge, 1996], we omit this term to maintain model simplicity, on the assumption that the lateral temperature gradient is negligible in comparison to the steep vertical temperature gradient. We assume a two-phase system. A rock phase, with void space equivalent to porosity, φ, is uplifted at a spatially constant vertical rock advection velocity, v u , relative to a fixed coordinate axis. A fluid phase saturates the void space and is uplifted with the rock, but also flows relative to the rock with vertical Darcy velocity (flux), q. These velocities combine to give an effective vertical advection velocity, v e . The combined fluid-rock medium has effective volumetric heat capacity, C e and effective thermal diffusivity, κ e [Buntebarth and Schopper, 1998;Eppelbaum et al., 2014]. Our models are based on the equation for steadystate heat flux, Q, in one dimension, z, representing depth below the surface: We assume a vertical profile comprising 1 to 4 sections within which q and Q are constant. Section boundaries represent intersections with structures that influence q in the adjacent sections, by modifying the magnitude or orientation of the 3D flow vector. Some such structures (e.g. fracture zones) may also transport fluid through the profile at very high local flow rates. Although we assume that lateral temperature gradients are small in comparison to vertical temperature gradients, along pathways of intense localized flow these gradients may be significant, such that through-flowing fluids act as an additional source or sink of heat. Material properties are fixed at constant values (see Supporting Information). v u is assumed constant at 8 mm·yr -1 [Little et al., 2005]. A best-fit solution (minimum root-mean-square deviation, RMSD) is calculated with q and Q as fitted parameters for each section. Model errors are 95% confidence intervals from model regression, including uncertainties in the fixed material properties and v u .
For initial comparison we consider two models for a single homogeneous section that simulate conditions of no fluid advection of heat. Model 0A simulates conduction only (v e ≈ 0), with Q set to a typical continental background value of 70 mW·m -2 [Stein, 1995;Davies and Davies, 2010]; Model 0B simulates rock advection and conductive heat exchange only (q = 0 and Q = 358 mW·m -2 based on uplift of rock and water at 500 °C above surface temperature from 15 km depth [Allis et al. 1979]. Both models poorly represent the measured data (Fig. 2a, RMSD = 27.47 °C and 11.01 °C respectively) indicating that both fluid flow and rock advection are important. In comparison, our Model 1, comprising one homogeneous section with constant q and Q, gives a visibly good fit to data (RMSD = 1.67 °C), although it diverges at the three distinct discontinuities in gradient (Fig. 2b).
To investigate these divergences the model domain is split into four sections in Models 2 and 3. Section boundaries are positioned at the three major discontinuities in temperature gradient. As such model sections (A1, A2, B, C) correspond to Layers A (deformed Alpine Schist), B (lacustrine sediments) and C (alluvial gravels), with Layer A further split into A1 and A2 sections at 691 m depth, based on the change in geothermal gradient at this depth. Model 2 allows q to vary between sections, with Q in each section dependent on q. This model does not allow for modification of heat flux by through-flowing fluid at section boundaries. Model 2 gives a strong fit (RMSD = 0.66 °C) and also reproduces the discontinuities in gradient that Model 1 cannot. In Model 3, Q is allowed to vary between sections, independent of q. As such this model allows for additional heat sources or sinks associated with intense fluid through-flow along any highly permeable pathways at section boundaries. Model 3, with RMSD = 0.50 °C, produces a tighter fit to discontinuities than Model 2, but at the expense of increased model complexity with three additional fitted parameters. Figure 2. a. DFDP-2B stratigraphy, borehole head and temperature data [Sutherland et al., 2017], with the results of single-section models based on typical continental conductive heat flow (Model 0A) and on heat transport by 1D rock advection (Model 0B). b. Configurations and results for single-section Model 1 and four-section Models 2-3. For Model 2 Q changes in proportion to q. For Model 3 Q is allowed to vary independently from q. In all models conductive heat flux gradually increases and advective heat flux decreases with decreasing depth, as fluid and rock approach surface temperature.
Step changes in heat flux associated with q changes at section boundaries are superimposed upon this trend in multi-section models.

Results
We first consider generalized, large-scale insights into heat transport near the Alpine Fault. The poor match of models neglecting fluid flow (Models 0A and 0B) to the DFDP-2B thermal profile highlights the importance of fluid heat transport (Fig 2). Although Models 0A and 0B could provide a better fit if Q was set to a much higher value, there is no obvious physical basis for this in the context of only considering conduction and rock advection.
Best-fit upward Darcy velocities for Models 1-3 range from 5.3×10 -10 m·s -1 to 2.8×10 -8 m·s -1 (Fig. 2b). Total model heat flux across the observed interval is ~720 mW·m -2 (Model 1). This is higher than published values for other boreholes in South Island [Shi et al., 1996], with the closest heat flux a value of 320 mW·m -2 measured in the Waiho Valley, near Franz Josef. In addition to exceeding the heat flux that can be produced by simplified onedimensional rock advection (Model 0B), this Model 1 heat flux value is also significantly greater than values (<300 mW·m -2 ) predicted by more sophisticated rock advection models [Shi et al., 1996;Upton et al., 2011]. The heat flux associated with fluid flow is more than three times larger than that associated with rock advection, for all model best-fit values (Fig.  2b). The Péclet number for the modelled interval, based on heat advection by fluid flow alone is: where is the model section length. Model 1 has of ~1, indicating that heat transport by fluid flow is of comparable significance to conduction over the length of the borehole.
The visibly good fit of Model 1 indicates that the assumption of constant q across the drilled interval is suitable for thermal modelling at this scale. Models 2 and 3 improve this fit, by allowing for variations in q, which show no trend with depth and are thus likely to result from localized features. Models 2 and 3 offer insights into how fluid and heat flow is modified by local geological structures. Section A1 corresponds to the deeper part of the metamorphic basement, where thermal gradient is shallowest (~31 °C·km -1 between ends of the section). The geothermal gradient increases to ~139 °C·km -1 in section A2. The best-fit Model 2 infers a fluid sink between these sections to explain this change in gradient (Fig. 2b). Although the best-fit Model 3 predicts fluid inflow with a heat source at this boundary, solutions involving outflow also lie within the 95 % confidence interval of this model and have RSME equal to the best-fit solution to two decimal places (see Supporting Information). A highly fractured zone, intersecting DFDP-2B at this depth, could laterally remove up-flowing fluids from the profile. It could also act as a heat source, transporting hotter fluids to the model profile. The observation of a radon spike at the A1-A2 boundary depth during drilling is consistent with the penetration of a highly permeable pathway such as a strongly fractured zone at this depth. Alternatively, intersection with a barrier to flow, such as the slip plane of a minor fault could also divert up-flowing fluids laterally. Given the sharpness of this gradient discontinuity and the orientation of the borehole approximately normal to mylonitic foliation at this depth, the fractured zone or fault is likely to be perpendicular to the borehole and therefore near parallel to foliation and the Alpine Fault.
The geothermal gradient in section B, ~184 °C·km -1 , is higher than in section A2. Best-fit values from Models 2 and 3 give slightly higher upward q of 1.5×10 -9 m·s -1 to 4.0×10 -9 m·s -1 in section B compared 5.3×10 -10 m·s -1 to 6.2×10 -10 m·s -1 in section A2. Applying Darcy's Law to hydraulic data (Fig. 2a) and Darcy velocities from Models 2 and 3 give hydraulic conductivity, K ≥ 10 -8 m·s -1 for section B and K = 10 -9 m·s -1 for section A2. For water properties based on the section midpoint temperatures [Lemmon et al., 2017], this represents permeability, k ≥ 10 -15 m 2 for sediments in section B and k = 10 -16 m 2 for basement rocks in section A2. The increased q in section B could result from an influx of fluid from a profile elsewhere in the valley, where up-flow rates in the basement are higher. In this interpretation, fluids reaching the more permeable sediments flow laterally, down a hydraulic gradient to the observed profile. Lateral in-flow of fluid on shallow flow paths into the sediments at the valley flanks could also increase q. Model 3 results infer an additional heat sink, as might result from concentrated flow of cool water near the sediment-basement interface from the valley flank toward its centre.
Section C has a geothermal gradient of ~73 °C·km -1 , significantly lower than section B. Model 2 and 3 infer downward fluid and associated heat flow in section C to fit this. The alluvial gravels of section C are likely to have a high K, 10 -7 m·s -1 to 10 -2 m·s -1 [Domenico and Schwartz, 1990] and downward fluid fluxes may result from locally recharged water flowing toward the valley centre, potentially with 3D Darcy velocity vectors which are less steeply inclined than in sections below. At the boundary between Sections B and C, Models 2 and 3 predict convergence of down-flowing and up-flow fluids, both of which must travel laterally before reaching the surface. The best-fit Model 3 solution infers an additional heat sink at this boundary, suggesting that any focused lateral flow along this boundary occurs up a temperature gradient (e.g., towards the valley centre).

Discussion
1D models of the DFDP-2B temperature profile estimate a total heat flux of >700 mW•m -2 . This greatly exceeds that which can be generated by rock advection alone and indicates vertical Darcy velocities approximating to 7.8×10 -10 m·s -1 where parameters are assumed homogeneous, and potentially varying across geological structures in the range 0 to 8.4×10 -8 m·s -1 . Some previous models for heat flow in the Southern Alps have discounted the role of groundwater flow in heat transport, based on the small heat contributions associated with the few known warm springs that are rare localized features within a regional groundwater circulation [Allis and Shi, 1995]. However, most basement fluid upwelling is obscured by thick recent sediment fills of the valleys, and inflections in the temperature profile indicate fluid flow and heat advection around the basement-sediment boundary.
Single-section Model 1 indicates that, considering a depth averaged best fit q, heat transport by fluid flow is more than three times that by rock advection. Considering the full uncertainty range for q in all sections across Models 1-3, heat transport by fluid flow may vary with depth from >500 times greater than the heat transport by rock advection, to as low as zero, given that q = 0 m·s -1 lies within the uncertainty interval for Section A2 in Model 3, which includes both up-flow and down-flow. A Péclet number of ~1 for advection by fluid flow indicates it is also of comparable significance to conduction across the observed interval. Thus, heat advection by fluid flow, rather than being negligible, is likely to act as a major control on temperature distribution in the shallow crust, where high Darcy velocities persist and is key to understanding the elevated heat flux at DFDP-2B.
The Alpine Fault extends to ~35 km [Stern et al., 2007], with rock being advected from deep in the ductile crust. Reduction in fluid fluxes with depth [Sims et al., 2015;Menzies et al., 2016] means that although heat advected by fluid flow is much greater than that advected by rock advection at DFDP-2B, at some point much shallower than 35 km, heat transfer by fluid flow is likely to become negligible. Thus, rock advection supplies heat to the base of a shallow groundwater circulation system, above which fluid flow rates are high enough to advect a greater amount of heat than rock advection (Fig. 3). In this groundwater system, heat is transported from ridges to valleys. Water up-flow increases the upward heat flux beneath valleys, in excess of that supplied by rock advection and water down-flow decreases it beneath ridges. Heat flow data for the central Southern Alps are sparse, with a borehole in the Waiho valley [Shi et al., 1996] and the DFDP-1B [Sutherland et al., 2012] and 2B [Sutherland et al., 2017] boreholes providing the only published data. This data is consistent with the conceptual model proposed here, showing a significantly lower temperature gradient at DFDP-1B (63°C•km -1 ), which is not located in a major valley, compared to the other boreholes (≥95°C•km -1 ), which are. Thus, the high heat flux at DFDP-2B indicated by models depends upon both fluid flow and rock advection. Our results corroborate the finding of Sutherland et al. [2017], that the high geothermal gradient at DFDP-2B can be attributed to heat transport by both rock advection and fluid flow, while further indicating that the latter is the more significant advector of heat in a shallow groundwater system beneath the Southern Alps extending at least to the base of DFDP-2B. Discontinuities in the measured temperature gradient can be explained simply in terms of changes in q. Changes in thermal properties at lithologic boundaries may also contribute to discontinuities, but are not required to reproduce them.
Although there is no evidence of a systematic decrease in q across the observed interval in Models 2 and 3, extrapolation of high Darcy velocities at DFDP-2B is limited at the deepest to the Alpine Fault (≤1.5 km below the surface at this site). Beneath upstream sections of valleys, where the Alpine Fault is deeper, high upward Darcy velocities could persist deeper than at DFDP-2B. Fluid flow will remain the dominant mode of vertical heat transfer where q is >2×10 -10 m·s -1 . Vein isotopic data indicates that Darcy velocities of this order may persist as deep as the brittle-ductile transition at ~6 km [Menzies et al., 2016].

Conclusions
For the simplest model (single-section Model 1), an upwards Darcy velocity of 7.8×10 -10 m·s -1 provides a good approximation to measured temperatures at DFDP-2B. Discontinuities in the measured temperature gradient can be understood in terms of variable fluid up-flow, resulting from flow into and out of the observed interval at major fracture zones and at the interfaces between sedimentary units. Modelled Darcy velocities, away from known warm springs, have magnitudes ≥5.3×10 -10 m·s -1 , which are sufficiently high that heat advection by groundwater is of comparable significance to conduction and more than three times greater than heat advection by rock advection across the DFDP-2B interval. The total heat flux at DFDP-2B is ~720 mW·m -2 . This is significantly greater than typical continental conductive heat fluxes of ~70 mW·m -2 [Stein, 1995;Davies and Davies, 2010] and also exceeds that predicted by models incorporating only conduction and rock advection (<300 mW·m -2 ). The high heat flux is interpreted as resulting from the regional groundwater flow system that is strongly influenced by the ridge-valley topography developed perpendicular to the Alpine Fault. This enhances the high upward heat flux from rapid rock advection beneath valleys, and reduces it beneath ridges. Thus, we infer that groundwater flow is a major control on the temperature field in the shallow (<6 km depth) Southern Alps. Groundwater-driven heat flow will also influence the temperature-dependent mechanical properties of the fault materials [Kohlstedt et al., 1995] of the Alpine Fault throughout most of the seismogenic crust in this region.