Using Thermal Springs to Quantify Deep Groundwater Flow and Its Thermal Footprint in the Alps and a Comparison With North American Orogens

The extent of deep groundwater flow in mountain belts and its thermal effects are uncertain. Here, we use a new database of discharge, temperature, and composition of thermal springs in the Alps to estimate the extent of deep groundwater flow and its contribution to the groundwater and heat budget. The results indicate that thermal springs are fed exclusively by meteoric water and make up 0.1% of the total groundwater budget. Spring water circulates on average to a depth of at least 2 km. The net heat extracted from the subsurface equals 1% of the background heat flow, which equals an average thermal footprint of 7 km2. Cooling by downward flow and heating by upward flow are three and two times higher than the net heat flow, respectively. Comparison with North American orogens shows that hydrothermal activity is higher in areas with high relief or areas under extension.

Groundwater flow can strongly influence temperatures in orogens (Bodri & Rybach, 1998;Forster & Smith, 1989;Maréchal et al., 1999). The extent of the thermal effects of groundwater flow has been a long-standing question in low-temperature thermochronology, where anomalous thermochronological ages are often caused by hydrothermal activity (Dempster & Persano, 2006;Duddy et al., 1994;Louis et al., 2019). Data from thermal springs have suggested that in parts of the Himalayas (Derry et al., 2009;Whipp & Ehlers, 2007) and the Cascades (Ingebritsen et al., 1992), groundwater may be responsible for more than 50% of the heat transport. Yet the overall contribution of groundwater flow to the heat budget of orogens has not been quantified systematically to our knowledge.
Thermal springs are the most visible surface expression of deep groundwater flow in orogens. Data on spring discharge and temperature offer opportunities to quantify deep groundwater flow and its thermal effects (Luijendijk, 2019;Manga, 2001). The only analysis of thermal springs in orogens known to us covers North America (Ferguson & Grasby, 2011). However, this study did not report the contribution of thermal springs to the groundwater and heat budget of orogens. Country-wide estimates of spring discharge and heat flow exist for Japan (Fukutomi, 1970;Muraoka et al., 2006;Sumi, 1980) but include a large number of pumping wells in addition to natural springs. Here, we report a new database of thermal springs in the Alps that includes temperature, discharge, hydrochemistry, and isotope data. We use this database to quantify groundwater origin, circulation depth, and the contribution of thermal springs to the groundwater and heat budget of the Alps. In addition, we compare the results to published data for orogens in North America and explore the relation between hydrothermal activity and relief, groundwater recharge, and background heat flow.

Methods
Our database consists of 311 thermal springs in the Alps, which includes 210 springs with discharge data, 94 springs with hydrochemistry data, and 92 springs with stable isotope ( 18 O and 2 H) data, as described in Supporting Information S1. The data set predominantly consists of natural springs but also includes 23 springs that are tapped by wells that are less than 100 m deep. The data set is available on Pangaea (Luijendijk et al., 2020a). We consider springs thermal springs if their temperature is at least 5°C higher than the mean annual temperature of the area where the spring is located (White, 1957). We used the watershed that each spring is located in to define the area.
The origin of spring water was assessed by comparing the δ 18 O and δ 2 H values in thermal springs to precipitation (IAEA & WMO, 2019) and values for water derived from magmatic devolatilization and metamorphic dehydration (Sheppard, 1986;Yardley, 2009). Groundwater circulation temperature was estimated using a compilation of empirical quartz geothermometer equations (Verma et al., 2008) ( Figure S1). Circulation depth was calculated by dividing the circulation temperature with the average geothermal gradient, following Grasby and Hutcheon (2001). Note that this assumes that groundwater flow does not affect the geothermal gradient strongly, which may not be true in some cases. The calculated circulation temperatures and depths are minimum estimates, because empirical geothermometer equations may underestimate circulation temperatures as they do not take into account mixing effects and precipitation of SiO 2 during upward flow (Ferguson et al., 2009). Mixing between shallow and deep sources has been demonstrated for a number of springs in the Alps and may be relatively common (Buser et al., 2013;Sonney et al., 2012;Waber et al., 2017). Maximum circulation temperatures of the deep groundwater component may be double that of bulk samples (Diamond et al., 2018). The size of the contributing area for each spring was calculated by dividing spring discharge by recharge rate at the spring location (De Graaf et al., 2015).
Heat exchange between the subsurface and the groundwater that discharges in springs, as illustrated in Figure S2, is given by where H is the heat gain or loss of water (W), ρ f is the water density (kg m −3 ), c f is specific heat capacity of water (J kg −1 K −1 ), T 1 and T 2 are groundwater temperature at the start and end of the flow path (K), and Q is groundwater discharge (m 3 s −1 ). The term H v represents the heating of water caused by viscous dissipation, which can be important in relatively cold springs with a high relief in their contributing area (Manga & Kirchner, 2004). Viscous dissipation is given by where h 1 and h 2 are the hydraulic head at the start and the end of the flow path (m). Note that for flow between the recharge point and the spring, Equation 2 is equivalent to the equation by Manga and Kirchner (2004) that includes the elevation difference instead of the difference in hydraulic head (h 2 − h 1 ). However, for flow to and from circulation depth one needs to take into account the conversion of gravitational potential to and from elastic potential to account for the change in fluid pressure. We therefore use hydraulic head in Equation 2, which represents the potential energy of groundwater per unit weight (Hubbert, 1940). Note that this equation is a simplification that assumes that density is independent of fluid temperature. This introduces an error that is relatively small given the high uncertainty of the other terms in this equation.
The net heat exchange with the subsurface is calculated by substituting T 1 by the average recharge temperature T r (K) and by substituting T 2 by the spring temperature T s (K) in Equation 1. For the viscous dissipation term, h 1 is substituted by the average elevation of the recharge area z r (m), and h 2 is substituted by the elevation of the spring, z s (m). For heat exchange during downward flow, we use the average recharge temperature T r (K) as a value for T 1 and the circulation temperature T m (K) as the value for T 2 in Equation 1. For the heat exchanged during upward flow, T m is used for T 1 , and the spring temperature T s is used for T 2. For the calculation of the viscous dissipation using Equation 2, h 1 is substituted by z r , and h 2 is substituted by h m , the hydraulic head at the circulation depth. For the calculation of viscous dissipation during upward flow h 1 is substituted by h m and h 2 by z s .
The values of average recharge elevation (z r ) and temperature (T r ) are unknown. However, the bounds of these parameters can be estimated. The maximum elevation was equal to the highest elevation in the watershed that the spring was located in, as calculated using watershed (Lehner & Grill, 2013) and elevation data (Danielson & Gesch, 2011). The minimum estimate for recharge elevation was equal to the elevation of the spring. Similarly, the minimum and maximum bounds of the recharge temperature were calculated as the recharge temperature at the spring location and the recharge temperature at the highest elevation in the watershed. Recharge temperature was assumed to equal average annual air temperature (Fick & Hijmans, 2017). The hydraulic head at the circulation depth (h m ) is unknown and was chosen as the mean of the spring elevation and the maximum recharge elevation. Note that the uncertainty in h m only affects viscous dissipation, which is a very small term for most springs.
We compared the heat flux of thermal springs with published background conductive heat flow data in the Alps (Global Heat Flow Compilation Group, 2013), with a mean background heat flow density of 76 mW m −2 and a mean geothermal gradient of 26°C km −1 (Text S2 and Figure S3). These values agree with a recently published thermal model of the Alps (Spooner et al., 2020). Given the low amount of data and the low variation of heat flow (Spooner et al., 2020), we use a single value for heat flow and geothermal gradient. The thermal footprint of springs was calculated by dividing the net heat flux by the mean background heat flow density.
The net heat flux of springs in the Alps was compared with several mountain belts in North America using a published database of spring heat flow (Ferguson & Grasby, 2011). These values did not take into account viscous dissipation, and heat flow may therefore have been slightly overestimated. In addition, we compared the spring discharge and heat flux with average groundwater recharge (De Graaf et al., 2015), relief, and background heat flow of each orogen. We used the standard deviation of elevation (Danielson & Gesch, 2011) as an indicator for relief. The background heat flow for North America was obtained using natural neighbor interpolation (Sibson, 1981) of heat flow data (Blackwell et al., 2011;Blackwell & Richards, 2004) ( Figure S4).
Discharge, SiO 2 and 2 H/ 18 O data were not available for all of the springs in the data set (Table S1). We explored whether this could introduce bias in our analysis by comparing temperature and discharge for the entire data set and for springs without discharge, SiO 2 , and 2 H/ 18 O data ( Figure S5) and found no significant difference between the distributions, which suggests that the distribution of data is not biased toward springs with low or high temperature or discharge.

Results and Discussion
The thermal springs are relatively evenly distributed over the Alps, with an aerial density of one spring per 650 km 2 and an average distance between springs of 14 km. The mean temperature of the springs is 22°C. The warmest springs are located in northern Italy (Figure 1a), where the La Bollente spring 10.1029/2020GL090134 Geophysical Research Letters (Bortolami et al., 1983) and Fonte Boiola (Waring & Blankenship, 1965) measured 71°C and 70°C, respectively. The highest rates of water discharge were observed in the eastern Alps (Figure 1c), where the Fontanon spring in Italy (Cantonati & Spitale, 2009) and the Kroparica spring in Slovenia (Philipp, 2015) discharge 0.22 and 0.31 m 3 s −1 , respectively. The locations, temperatures, and discharge magnitudes of thermal springs in the Alps are shown in Figure 1.
The mean discharge rate of 210 thermal springs with discharge data is 8.9 × 10 −3 m 3 s −1 , and total discharge is 1.9 m 3 s −1 , which equals to an annual rate of 0.3 mm a −1 when spread out evenly over the Alps. This total spring flow is only 0.06% of the average recharge rate in the Alps of 500 mm a −1 (De Graaf et al., 2015). This means that a large majority of groundwater in the Alps follows shallow flow paths and discharges in rivers or lakes instead of in springs.
The spring water for 92 springs with stable isotope data plot within 1‰ δ 18 O value of the global meteoric water line and local precipitation values, which indicates a meteoric origin of thermal spring waters (Figure 2a). The compiled SiO 2 concentrations and the calculated values of circulation temperature and depth for 91 of the thermal springs are shown in Figure 2. The circulation temperatures range up to 130°C , which corresponds to a depth of 4,700 m below the surface based on average conductive geothermal gradients in the Alps (Text S2; Figure S3). The mean circulation temperature is 65°C, which is equal to a depth of 2,300 m below the surface. Note that these depths are minimum estimates due to the limitations of the SiO 2 geothermometer discussed in section 2.
The calculated net heat flux of these mapped thermal springs in the Alps is shown in Figure 3. The average net heat flux for thermal springs equals 0.4 to 0.7 MW. The reported range of heat flux and thermal footprint represents end-member (minimum and maximum) uncertainty estimates. Note that the net heat flux as defined here is equal to heat extracted from the subsurface and does not include viscous dissipation, which is discussed more detail below. The two springs with the highest net heat flux are the Urquelle in Villach, Austria, with a net heat flux of 7.0 to 11 MW and Lavey-les-Bains in Switzerland with a net heat flux of 5.3 to 7.5 MW. The calculated net heat flux for 210 thermal springs for which the average temperature and discharge rate are known equals 84 to 146 MW, which equals a spatially averaged heat flow density for the entire Alps of 0.42 to 0.72 mW m −2 . The heat flux by thermal springs is equal to 0.5% to 0.9% of the background heat flow. The thermal footprint of the springs is on average 5.3 to 9.1 km 2 (Figure 3d). This is much larger than the contributing area for each spring, which is on average 0.6 km 2 (Figure 3b).
Viscous dissipation contributes 0% to 26% of the total heat output of thermal springs in the Alps. The high uncertainty of the contribution of viscous dissipation is due to the unknown elevation and hydraulic head in the contributing area of springs as discussed in section 2. Viscous dissipation is the highest in low-temperature springs and is likely the dominant heating process for most springs with temperatures lower than 15°C ( Figure S6).
The net heat flux of the thermal springs is a combination of the cooling of the subsurface by downward groundwater flow and heating by upward groundwater flow ( Figure S2). Cooling by downward  (Yardley, 2009). VSMOW denotes Vienna Standard Mean Ocean Water, GMWL denotes the global meteoric water line (Craig, 1961), (b) SiO 2 concentration in thermal springs, (c) calculated circulation temperature, and (d) circulation depth. Circulation depth was calculated using an average geothermal gradient of 26°C km −1 . The min., max., and best values in panels c and d represent the circulation temperature and depth using the lowest, highest, and mean value of several geothermometer equations, respectively, as explained in section 2. Note that for three springs the SiO 2 value was lower than 5 mg L −1 , which is the lower limit for which the geothermometer equations are valid.

Geophysical Research Letters
LUIJENDIJK ET AL.
groundwater flow for 58 springs with data on circulation temperature is equal to 51 to 87 MW, and heating by upward flow for these springs is calculated as 30 to 56 MW. In comparison, for these springs the net heat flux (difference between the heat flux by downward and upward flow) is estimated as 20 to 32 MW, which illustrates the potential for advection to redistribute heat and impact subsurface temperatures, even when the net heat flux (i.e., the heat flux observed at the discharge point) is relatively small.
Comparison with orogens in North America shows that the thermal spring discharge and the spring heat flux of the Alps are, relatively, high (Figure 4). Note that in contrast to the Alps, the spring heat flux of other orogens may be overestimated because viscous dissipation was not taken into account. However, the differences in hydrothermal heat flux between the Alps and other orogens are higher than the relatively small contribution of viscous dissipation. We hypothesize that the differences in hydrothermal activity between the orogens can be explained by several factors. First, the Alps have a relatively high spring density (Figure 4c). Due to a lack of systematic exploration and mapping of springs, it is unknown if this represents a real difference in hydrothermal activity or whether this simply represents a more complete database for the Alps. Second, the hydrothermal activity appears to be associated with topographic relief (Figure 4d). Relief controls the maximum hydraulic gradient and is a measure of the driving force available for deep groundwater flow (Forster & Smith, 1988). Third, in orogens cut by fault systems presently in transtension and extension, fault permeability tends to promote deep groundwater circulation (De Paola et al., 2007;Grasby & Hutcheon, 2001;Hartman et al., 2018;von Hagke et al., 2019). The Sierra Nevada orogen shows the highest hydrothermal activity (Figures 4d-4f). Apart from the relatively high relief, this is likely caused by fault systems in the Sierra Nevada being dominated by active transtension (Bennett, 2011). The comparison between background heat flow and hydrothermal activity (Figure 4f) confirms earlier work that showed

10.1029/2020GL090134
Geophysical Research Letters that background heat flow does not control spring heat flux (Ferguson & Grasby, 2011). Note, however, that due to poor data coverage, heat flow in the McKenzie mountains, Coastal Mountains, and part of the Rocky mountains is highly uncertain ( Figure S4). Instead, hydrothermal activity in the western United States (Figure 4a) is most likely promoted by extensional tectonics and high strain rates that provide deep permeable flow paths along fault zones (Faulds et al., 2012). Note that in addition to the factors explored here, we anticipate that other factors such as lithology, deformation history, fault density, and structural style may play a role in explaining thermal spring activity in orogens.
The contribution of thermal springs reported here is a minimum estimate for deep groundwater flow and its thermal effects. The net heat flux by thermal springs is probably higher because 32% of the thermal springs in the database had no discharge data. If one assumes that the distribution of temperature and discharge in the remaining springs would be the same as for the springs with data, then the net heat flux would be higher and equal approximately 1% to 2% of the background heat flow. In addition, thermal springs likely represent only a part of the total deep groundwater flow, because a large part of deep groundwater may form blind hydrothermal systems without associated thermal springs at the land surface. Large quantities of groundwater usually discharge diffusely into rivers or lakes and may contain waters of hydrothermal origin, but relatively few estimates for the rate of spring versus diffuse discharge have been published. Fluid and heat budget calculations in a river basin in the Himalayas (Derry et al., 2009), the Beowawe system (Olmsted & Rush, 1987), and the Borax lake geothermal area (Fairley & Hinds, 2004) in the Basin and Range Province showed that

10.1029/2020GL090134
Geophysical Research Letters thermal spring discharge equals at most 1/3 of the total hydrothermal discharge. In the Alps, there is subsurface in situ evidence for diffuse thermal water discharge in tunnels and boreholes. For example, the construction of a 75 m-long tunnel in Font Salée resulted in discharge with an estimated heat output of 69 MW (Silvestre, 1991) that exceeds any of the natural springs in the Alps. Similarly, in Brigerbad (Buser et al., 2013) and La Léchère (Thiébaud et al., 2010), shallow (<100 m) boreholes have tapped artesian groundwater with a higher heat output than existing natural springs at the same locations.

Summary and Conclusions
We present and analyze a new database of thermal springs in the Alps. Our analyses suggest that thermal springs are fed by meteoric water and the discharge of the springs comprises 0.1% of the total meteoric groundwater budget of the Alps. The circulation depth of the water that discharges in the thermal springs ranges up to 4,700 m and has a mean value of 2,300 m. This is an underestimate due to the limitations of the SiO 2 geothermometer on which these numbers are based (Ferguson et al., 2009), which means that meteoric groundwater circulates in a substantial part of the upper crust. The net heat flux by thermal springs is approximately 0.5% to 0.9% of the background heat flow of the Alps and is a combination of two larger terms: cooling by downward groundwater flow and heating by upward groundwater flow, for which the magnitude is approximately three and two times the net heat flux, respectively. These terms are probably underestimated due to the underestimation of circulation temperatures. Spring heat flux is also underestimated because no discharge data were available for 32% of the thermal springs in the Alps. In addition, given that not all deep groundwater flow discharges in thermal springs, the estimates reported here are a minimum estimate for total amount of deep groundwater flow and its contribution to the heat budget.
Although the contribution to the total heat budget is small, hydrothermal activity can locally be important for subsurface temperatures and heat flow. The average thermal footprint of springs is 7 km 2 , which means that on average heat flow and temperatures are affected in an area with a radius of 1500 m around springs, and larger if one assumes that springs do not harvest all the background heat flow over a particular area but harvest a part of the background heat flow over a larger area. Comparison with orogens in North America shows that hydrothermal activity is highest in orogens with high relief and areas undergoing extension or transtension.

Data Availability Statement
The thermal spring database has been published at Pangaea (Luijendijk et al., 2020a). Note that this database contains two data sets, one that includes deep wells and areas close to the Alps (Luijendijk et al., 2020b) and a second file that includes only the data used in this study (Luijendijk et al., 2020c). The notebooks used for data analysis have been published at Zenodo (Luijendijk, 2020) and are also available on GitHub (https:// github.com/ElcoLuijendijk/thermal_springs_alps).