Beo v1.0: Numerical model of heat flow and low-temperature thermochronology in hydrothermal systems

Low-temperature thermochronology can provide records of the thermal history of the upper crust and can be a valuable tool to quantify the history of hydrothermal systems. However, existing model codes of heat flow around hydrothermal systems do not include low-temperature thermochronometer age predictions. Here I present a new model code that simulates thermal history around hydrothermal systems on geological timescales. The modelled thermal histories are used to calculate apatite (U-Th)/He (AHe) ages, which is a thermochronometer that is sensitive to temperatures up to 70 °C. The modelled AHe 5 ages can be compared to measured values in surface outcrops or borehole samples to quantify the history of hydrothermal activity. Heat flux at the land surface is based on equations of latent and sensible heat flux, which allows more realistic land surface and spring temperatures than models that use simplified boundary conditions. Instead of simulating fully coupled fluid and heat flow, the code only simulates advective and conductive heat flow, with the rate of advective fluid flux specified by the user. This relatively simple setup is computationally efficient and allows running larger numbers of models to quantify model 10 sensitivity and uncertainty. Example case studies demonstrate the sensitivity of hot spring temperatures to the depth, width and angle of permeable fault zones, and the effect of hydrothermal activity on AHe ages in surface outcrops and at depth.

. While numerous model studies have addressed fluid and heat flow in terrestrial hydrothermal systems on geological timescales (Wieck et al., 1995;Banerjee et al., 2011;Howald et al., 2015;Volpi et al., 2017), the effect of hydrothermal activity on thermochronometer ages has to my knowledge not been modelled. The sole exceptions are Person et al. (2008), who found good agreement between apatite fission track ages around the Carlin gold deposit and ages predicted by numerical models, and Luijendijk (2012), who combined an advective and conductive heat flow model, which 5 was a precursor to the model code presented here, to model heat flow and apatite fission track ages around a hydrothermally active normal fault.
While coupled fluid and heat flow models can provide realistic reconstructions of the thermal history of hydrothermal systems, they are also relatively computationally expensive, which may limit the possibility to explore the response of these systems to different parameters. Heat flow data or thermochronology data are frequently scarce and relatively uncertain. These 10 data can often be explained by a number of different parameter combinations on for instance the age, duration and flow rates of hydrothermal systems. Here I present a new advective-conductive heat flow code, Beo, that can be used to model heat flow and apatite (U-Th)/He thermochronometer ages around hydrothermal systems over geological timescales. In contrast to coupled fluid and heat flow models, fluid flux is prescribed and is therefore not a function of the driving forces of fluid flow and permeability of the subsurface. This makes the code relatively computationally efficient and enables running larger numbers 15 of model experiments with variable parameters. In contrast to inverse thermal models such as HeFTy (Ketcham et al., 2007) and QTQt (Gallagher, 2012) which reconstruct thermal histories of one or more samples, Beo models 2D temperature fields over time, which can be compared to multiple low-temperature thermochronology samples and which allows testing different hypotheses for the history and characteristics of hydrothermal activity.
2 Model development 20 Beo was designed to model advective heat flow in and around a single main fluid conduit. A typical models setup is shown in Fig. 1. The model domain contains a main fluid conduit, which represents permeable fault zone. The main conduit is attached to one or more horizontal fluid conduits that can be used to model lateral flow in and out of permeable formations, that for instance represent alluvial sediments or permeable fractured rocks.
The following sections describe the equations used by Beo to model subsurface and surface heat flux and the apatite (U- 25 Th)/He thermochronometer in hydrothermal systems.

Advective and conductive heat flow
The heat flow equation used by Beo to model conductive and advective heat flow in the subsurface is given by: The remaining fluid discharges at the surface. Heat transfer between the surface and the atmosphere is modelled as a conductive heat flow, with a variable thermal conductivity based on equations for sensible and latent heat flux. Heat transfer in the subsurface is determined by the specific thermal conductivities of the local lithologies. The land surface is lowered over time to account for erosion.
, and q is fluid flux (ms −1 ). Subscripts b and f denote properties of the bulk material and the pore fluid, respectively. Beo solves the implicit form of the heat flow equation by discretization of the derivative of T over time: where ∆t is the size of a single timestep (s), T t (C) denotes the known current temperature and T t+1 (C) denotes the heat flow equation, i.e. the heat flow equation with the term ∂T ∂t and the heat advection term (q) set to zero. Beo uses the generic finite element code Escript (Gross et al., 2007a(Gross et al., , b, 2008 to model heat transfer. Escript employs Python bindings to the internal c++ model code and enables the use of multiple processors to increase computational power. Mesh generation was performed using GMSH (Geuzaine and Remacle, 2009) using a Python interface included in Escript. The discretized heat flow equation was solved using the GMRES solver (Saad and Schultz, 1986).

Land surface heat flux
Typical thermal boundary conditions for subsurface heat and fluid flow models are either a specified temperature or a specified heat flux at the top model boundary, which is usually chosen as the land surface. However, in transient hydrothermal systems it is difficult to simulate realistic temperature using these boundary conditions, especially in cases where fluid discharges at the surface. Assigning a specified temperature or heat flux would require knowledge of the change in fluid temperatures over time, 10 which is only rarely available.
Temperatures at the land surface are predominantly determined by latent and sensible heat flux (Bateni and Entekhabi, 2012). Beo uses an approach that is to my knowledge new in hydrothermal model codes, and models the heat flux in a layer of air overlying the land surface. The top boundary of the model domain is located several meters in the air, and is assigned a specified temperature that reflects the average annual air temperature. Latent and sensible heat flux from the land surface are 15 approximated by assigning an artificially high value of thermal conductivity to the layer of air. The thermal conductivity of the air layer is calculated using equations for latent and sensible heat flux described below.
The motivation for including the series of equations below in the model code is to simulate realistic land surface and spring temperatures in transient models. Note that the implementation does not incur high computational demands, the equations are evaluated for land surface nodes only and in contrast to the heat flow equations in section 2.1 these equations can be solved 20 directly without numerical methods.
Following Bateni and Entekhabi (2012) the sensible heat flux at the land surface is given by: where H is the sensible heat flux (Wm −2 ) ρ is density (kg m −3 ), c a is the specific heat of air (J kg −1 K −1 ), r a is the aerodynamic resistance (s m −1 ), T a is the air temperature at a reference level (C) and T s is the surface temperature (C). This 25 expression can be combined with Fourier's law: where q is heat flux (Wm −2 ), K is thermal conductivity (Wm −1 K −1 ), ∆T = (T a − T s ) (K). Combining equation 3 and 4 and equalling H and q yields an expression for the effective thermal conductivity (K s ) between the surface and the reference level z: where ∆z is the difference between the surface and the reference level (m).

5
Latent heat flux is given by (Bateni and Entekhabi, 2012): where LE is the latent heat flux (Wm −2 ), ρ a is the density of air (kg m −3 ), L is the specific latent heat of vaporization (J kg −1 ), which is 334000 Jkg − 1, q s is the saturated specific humidity at the surface temperature kg kg −1 ), q a is the humidity of the air (kg kg −1 ). Combining this with Fourier's law gives the heat transfer coefficient for latent heat flux (K l ) as: The saturated specific humidity (q s ) was calculated as (Monteith, 1981): where e s is saturated air vapour pressure (Pa), P a is surface air pressure (Pa). The saturated air vapour pressure was calculated using the Magnus equation (Alduchov and Eskridge, 1996): Air pressure was assumed to be 1 × 10 5 Pa. The thermal conductivity assigned in the air layer is the sum of the heat transfer coefficient for latent heat flux (K l ) and sensible heat flux (K s ).
The resulting heat flux at the land surface is predominantly a function of the aerodynamic resistance (r a ). We use a value of 80 s m −1 following values reported for areas covered by short vegetation (Liu et al., 2007), and use a range of ±30 s m −1 to 20 quantify the uncertainty of r a and its effect on modelled spring temperatures.
The calculated heat transfer coefficient shows a strong dependence on surface temperature and aerodynamic resistance, as shown in Fig. 2. This implies that transient numerical models of hydrothermal systems cannot use a fixed heat transfer coefficient at the land surface, because changes in land surface and spring temperature over time change the heat transfer coefficient.

Boiling temperature
The model code simulates conductive and advective heat flow for a single fluid phase. Beo contains an option to cap subsurface temperatures to the boiling temperature curve, which is calculated using a 3rd order polynomial fit to data by the National

Institute of Standards and Technology (2018):
T max = 3.866 log(P ) 3 + 25.151 log(P ) 2 + 103.28 log(P ) + 179.99 where T max is the maximum (boiling) temperature in the system (°C), and P is fluid pressure (Pa), which is assumed to be hydrostatic. Using the boiling temperature as upper limit ensures that subsurface temperatures remain realistic in a system where vapour is present. Note however that this approach is a simplification in which the latent heat of vaporization is ignored and in which bulk thermal conductivities are not adjusted for the presence of a vapour phase. Therefore, this is intended as a first order approximation of temperatures in multi-phase flow systems, but for more realistic models it would be preferable 10 to use multi-phase flow codes such as Hydrotherm (Hayba and Ingebritsen, 1994). Note that the choice for capping modelled temperature by the boiling temperature avoids the high computational expense of including phase transformations in the model code, which would necessitate solving an additional equation for each node and each timestep, and which would likely place more constraints on spatial discretization and timestep size. 15 For modelling systems that are active over longer timescales Beo can take into account erosion or sedimentation by lowering or raising the land surface over time. This is done in a stepwise fashion, with the default value set to steps of 1 m. The implementation of erosion is important for low-temperature thermochronometers, because it exposes rocks that have been buried deeper and may have experienced more hydrothermal heating than rocks at the surface that are buffered by surface temperatures.

Apatite (U-Th)/He thermochronology
The modelled temperature history was used to calculate the response in low-temperature thermochronometer apatite (U-Th)/He 5 (AHe). The AHe thermochronometer is based on the sensitivity of the diffusion of helium in apatite minerals to temperature.
At high temperatures helium diffuses out at a rate equal or higher than the production by radioactive decay, whereas at low temperatures helium is retained in the apatite minerals. The AHe thermochronometer is sensitive to temperatures ranging from approximately 40 to 70°C (Reiners et al., 2005). AHe ages were calculated by solving the helium production and diffusion equation for apatites using the Eigenmode method, following Meesters and Dunai (2002a) and Meesters and Dunai (2002b), 10 which is a computationally efficient method that provides the same results as more computationally demanding finite difference methods (Meesters and Dunai, 2002a, b). Helium production and diffusion is described by: where C is the concentration of helium (mol m −3 ), D is the diffusion coefficient of helium (m 2 s −1 )and U is the helium production rate (mol m −3 s −1 ). The term S p denotes the probability that an emitted alpha particle stops at location x,y,z in 15 the apatite crystal (dimensionless), which is required because the long stopping distance of alpha decay (approximately 21 µm) means that some of the helium that is produced by radioactive decay is not retained inside the crystal. The helium age is calculated using the modelled average concentration of helium in the crystal (C avg ) following: The helium diffusion model assumes a spherical shape of the apatite crystal, and can be compared to measured AHe ages 20 when using an equivalent spherical shape with the same surface to volume ratio as the actual crystal (Meesters and Dunai, 2002a).
The diffusivity (D) of helium in apatites depends on temperature and on radiation damage. Radiation damage slows helium diffusion (Shuster et al., 2006;Flowers et al., 2009;Gautheron et al., 2009). The effects of radiation damage on the diffusivity of helium is calculated using the RDAAM model by Flowers et al. (2009), in which the diffusivity is a function of the density 25 of damage tracks. The annealing of these tracks is dependent on temperature and is assumed to occur at the same rate as fission tracks following established models of fission track annealing (Ketcham et al., 2009(Ketcham et al., , 2007. Beo also includes an option to calculate helium diffusivity following equations by Farley (2000) or Wolf et al. (1998) instead that do not take radiation damage into account.
The parameter S p is used to correct He concentration in the apatite crystal for the chance that alpha particles that are ejected implementation of this parameter in the helium diffusion model. I used an alpha stopping distance of 21 µm (Ketcham et al., 2011).

Model verification
I first validated the transient conductive heat flow in the numerical model using an analytical solution for the cooling of an intrusive in the subsurface. The solution for temperature change of an initially perturbed temperature field is (Carslaw and 5 Jaeger, 1959): where T b is the background temperature (C), T i is the temperature of the intrusive (C), L is the length of the intrusive (m), x is distance from the intrusive (m) and κ is thermal diffusivity (m 2 s −1 ) Following Ehlers (2005) I use this equation to simulate cooling of an intrusive in the subsurface as a test case for the model 10 code. The numerical and analytical solutions for cooling are shown in Figure 3. The intrusive body has an initial temperature of 700°C, and stretches from 0 to 500 m distance. The background temperature is 50°C. The solutions match to within 1.0°C.
In addition, the performance of the model code was evaluated using an analytical solution of steady-state heat advection and conduction by Bredehoeft and Papaopulos (1965). The solution describes heat advection in a one-dimensional system with fixed temperatures at the top and bottom boundaries of the model domain. The analytical solution is given by: and where T is temperature (K), z is depth (m), ∆T is the temperature difference between the top and the bottom of the domain the domain (m) and K is thermal conductivity (W m −1 K −1 ). A comparison between the numerical solutions by Beo and the analytical solutions shows that the analytical and numerical solutions are identical ( Fig. 3 and 4).
4 Model usage

Model input
Beo uses a single python file that contains all input parameters. The input parameters are Python variables. The parameter file 25 is organized into two parts, which each part contained in a Python class, named ModelParams and ParameterRanges, respec-  combination, which can be used to explore parameter space.

Model output and visualization
Beo generates output to comma-separated files, model-specific output files using the Python pickle module, and output of the modelled temperatures and fluid fluxes in VTK format. The comma-separated files contain a copy of all input parameters for each model run, along with several statistics for the model output such as modelled average change in temperatures compared 10 to initial temperatures, modelled temperatures at the surface or user-specified depth slices, modelled AHe data and comparison to observed values. Modelled temperatures and fluxes can be saved as VTK files that can be used for model visualization using external software such as Paraview and Visit. In addition, the model results can be saved in a Beo specific file format that contains all modelled temperature and AHe data. These output files can be used by a separate script (make_figure.py) to automatically generate figures such the model results shown in this study (Fig. 5). 15

Application
The following section presents models of two active hydrothermal systems. The first case study consists of a model of the thermal history of the Baden and Schinznach spring system, a hydrothermal system and series of hot springs at the boundary of the Jura mountains and the Molasse Basin in Switzerland. This study demonstrates the potential and limitations of the use of spring temperatures and discharge to quantify the depth of fluid conduits, and the use of low-temperature thermochronology 20 to reconstruct the history of hydrothermal activity. The second case study is a model of borehole temperatures and borehole apatite (U-Th)/He data in a hydrothermal system at the Brigerbad spring in the Rhone Valley in the western Alps, which is based on data published by Valla et al. (2016). Note that the aim is not to provide detailed case studies, but to illustrate the possibilities of the model code. In addition to these examples, a separate study that uses the model code to quantify the history of the Beowawe hydrothermal system in the Basin and Range Province has been published as a preprint on EarthArxiv (Louis 25 et al., 2018) and is currently in press in Geology (Louis et al., in press]).

Model setup
The model of the Baden and Schinznach spring system is based on a model study by Griesser and Rybach (1989). The total heat flux in the Baden and Schinznach spring system that is used as a case study is 2.4 × 10 6 W, which was calculated using 30 spring discharge and temperature data reported by Sonney and Vuataz (2008) and an assumed recharge temperature of 10°C.
With a background heatflow of 0.07 Wm 2 the minimum contributing area for the heat output of the springs is 3.5 × 10 7 m 2 . This is slightly above the median value for springs in North America reported by Ferguson and Grasby (2011). The model study presented here therefore represents a terrestrial hydrothermal system with a relatively high, but not unusual heat output.
The area hosts a number of springs with a temperature of 30 to 47°C (Sonney and Vuataz, 2008), and an average discharge 5 along the fault of 2 × 10 −5 m 2 s −1 . Fluid flow is hosted in a relatively shallow thrust fault that dips around 50 degrees to a detachment level around 1000 m below the surface, which may be connected to a deeper normal fault (Griesser and Rybach, 1989;Malz et al., 2015).
The numerical model is based on a conceptual model shown in Fig. 1. The model only includes the discharge part of the hydrothermal system. Groundwater recharge is much more diffuse than discharge and has a negligible effect on subsurface 10 temperatures in comparison to focused groundwater discharge, as shown for instance by model experiments by Ferguson et al.

(2009).
A specified heat flow of 0.07 Wm −2 was chosen at the lower boundary (Griesser and Rybach, 1989). For the upper model boundary the air temperature is fixed at 10°C at an elevation of 2 m above the land surface, and the heat transfer at the land transfer is governed by sensible and latent heat flux following equations 3 to 9. Thermal conductivity was fixed at 2.5 The model results show a strong dependence of modelled spring temperatures on the assumed depth of the fluid conduit (Fig. 6a). The observed temperatures of the Baden and Schinznach spring system are only reached using fluid conduits that are 30 at least 7 km deep. In addition, comparison of a the effects of fluid conduits dipping 50 degrees and 65 degrees show a strong difference in spring temperatures (Fig. 6a). For a fluid conduit with a low dip angle the flow path from a particular depth is longer and therefore the amount of heat loss along the way is higher. Modelled spring temperatures for a fluid conduit of 50 degrees are too low to explain the observed temperatures in the system. The model results confirms the hypothesis proposed by Griesser and Rybach (1989) that the fluid source in the Baden and Schinznach hydrothermal system is likely a deep and steeply dipping normal fault that is connected to the more shallow thrust fault that hosts the springs near the surface.
In addition to the depth of the fluid conduit spring temperatures are sensitive to the assumed width of the fault zone (Fig. 6b).
The wider the fault zone, the lower the spring temperature. Note that in these model runs, the overall flux was kept at 2 × 10 −5 m 2 s −1 , which was redistributed evenly over the width of the fluid conduit. The narrower the fluid conduit, the higher the flow 5 velocity, and the lower the conductive heat loss along the way.
The model experiments also show a strong dependence of spring temperatures on the modelled heat flux at the land surface.
The key parameter governing latent and sensible heat flux at the land surface is the aerodynamic resistance (Fig. 2). The value of aerodynamic resistance strongly affects spring temperatures. Lower values of resistance, which correspond to more vegetated conditions (Liu et al., 2007), result in higher values of effective thermal conductivity and heat flux at the surface (Fig. 2), and 10 as a result lead to lower spring temperatures (Fig. 6b).

Hydrothermal activity and low-temperature thermochronology
The effect of hydrothermal activity on low-temperature thermochronology was explored by using modelled thermal history of the Baden and Schinznach hydrothermal system to calculate to calculate AHe ages. For models on longer timescales exhumation plays a role. Due to the buffering effect of air temperature rocks at deeper levels heat up much more than rocks 15 close to the land surface (Fig 5). The rate of exhumation therefore determines the strength of the hydrothermal perturbation of thermochronometer ages. The effect of exhumation rates is explored using two different exhumation rates, representing slowly exhuming areas such as passive margins (1 × 10 −4 m a −1 ) and moderate exhumation rates representing orogens ( . Note that to better compare the model results, the initial AHe age, i.e., the AHe age before the start of hydrothermal activity, was set to the same value for both model experiments. The results demonstrate that the effect of hydrothermal activity on thermochronometers is dependent on background ex-5 humation rates. For a model run with a high exhumation rate of 1 × 10 −3 m a −1 the width of the zone at the surface where samples are partially reset is 35 m after 15000 years, which is 10 m wider than the width of the fluid conduit (25 m) in this model experiment. This means that low-temperature thermochronometers can be used to quantify the history of active hydrothermal systems and hot springs, but only if sampling is very dense.
For samples located at 500 m depth thermochronometers are affected up to 850 m distance from the fluid conduit (Fig. 8).
10 This means that even over a relatively short timescale of one interglacial stage (∼15000 years), hydrothermal activity can affect low-temperature thermochronometers at the subsurface in large areas. The effect of hydrothermal activity may be important for the interpretation of thermochronometers from boreholes near hydrothermal systems or areas close to exhumed hydrothermal systems.

Model setup
The Brigerbad hydrothermal system contains three natural springs that are located at the valley floor at the transition from the crystalline bedrock of the Aar massif to low-permeable quaternary valley fill. The springs are fed by a mixture of shallow cold water and a deep groundwater source that is likely channelled upward along a major thrust fault that is located in the 5 valley floor and that forms the contact between the crystalline rocks of the Aar massif and low permeable sedimentary units of the Helvetic nappe (Buser et al., 2013;Valla et al., 2016). Subsurface temperatures in a 500 meter deep borehole close to the springs show a strongly elevated geothermal gradient of approximately 100°Ckm −1 that is caused by the upward flow of warm groundwater.
The hydrothermal system was modelled using a simplified setup that uses a wide fault zone that functions as a flow conduit 10 with a maximum depth of 3.5 km and an upper boundary at a depth of 100 m, where the fault is overlain by low-permeable Quaternary sediments (Buser et al., 2013). The depth of fluid flow is based on an estimated maximum fluid temperature of 110°C, which is based on geothermometers (Buser et al., 2013). The fault zone is attached to an inclined fluid conduit that channels fluids away laterally from the overlying Quaternary sediments towards the springs that are located at a distance of approximately 125 m north of the fault outcrop. In addition there is a second fluid conduit that channels shallow and cold 15 groundwater from the Aar massif to the north to the same springs and that was estimated to be 100 m deep following general permeability trends in fractured crystalline rocks (Ranjram et al., 2015;Achtziger-Zupančič et al., 2017). The total flux was estimated as 230 m 2 a −1 following reported discharge values of the springs (Sonney and Vuataz, 2008). Following published conceptual models and hydrochemical data (Buser et al., 2013)  The temperature history for each sample was calculated by converting exhumation to cooling rates using a geothermal gradient of 27°C km −1 , adding the initial background temperature based depth of the samples and subsequently combining the initial temperature history with the modelled temperature history, which covered the last 10,000 years of the thermal history of the samples.

Results
Initial model experiments that use a narrow flow conduit that corresponds to the typical width of a fault damage zone for a major fault of 100 m (Bense et al., 2013) fail to provide a good fit the the borehole temperature record and overpredict 5 temperatures at shallow depths (Fig. 9). The model only provides a good fit when deep groundwater is channelled upward along a relatively wide flow conduit, with a size of 500 m providing a good fit to the temperature data (Fig. 10). This confirms models put forward in earlier studies (Buser et al., 2013) that suggest broad distribution of permeable fractures and fluid flow, as opposed to flow that is confined to a narrow fault damage zone. The modelled temperatures do still show a misfit, which suggests that the flow paths may be even more distributed and complex than the simplified model shown here.

10
The thermal history results in modelled AHe ages that are largely unaffected by hydrothermal activity and that range from 1.4 Ma to 2.1 Ma (Fig. 10d). This is much higher than the measured values in the three samples reported by Valla et al. (2016). The modelled hydrothermal heating during 10,000 years is too short to affect AHe ages. The mismatch indicates that hydrothermal activity was active for much longer than the 10,000 years that were modelled here. The system is likely to have been active periodically, during interglacial stages.  match observed discharge rates and temperatures in systems that host thermal springs. The code can also quantify the thermal footprint of hydrothermal systems at the surface and at depth, and provides a tool to quantify the effects of hydrothermal activity on low-temperature thermochronometers. The effects of hydrothermal activity on thermochronometers depends strongly on the duration of activity, and the code provides new opportunities to use thermochronometers to quantify the history of hydrothermal systems and hot springs.

10
Code availability. The source code of Beo version 1.0 has been published at Zenodo (Luijendijk, 2018) and is accessible online (https: //zenodo.org/record/2527845). The source code is also available at a GitHub repository (https://github.com/ElcoLuijendijk/beo). The code is distributed under the GNU General Public License, version 3. The repository contains a readme file with a brief description of model installation, usage and output and a more extensive manual that includes a detailed description of all variables in the input data files. In addition, the repository contains several Jupyter notebooks that can be used to reproduce the model benchmarks discussed in section 3, and the input files for the model examples in section 5. Beo v1.0 depends on the generic finite element code escript (https://launchpad.net/ escript-finley) and the Python modules numpy, scipy and pandas. Braun, J., van der Beek, P., Valla, P., Robert, X., Herman, F., Glotzbach, C., Pedersen, V., Perry, C., Simon-Labric, T., and Prigent, C.: Quantifying rates of landscape evolution and tectonic processes by thermochronology and numerical modeling of crustal heat transport