Actual measurements of advective heat flux from Yellowstone hot springs (Wyoming, USA) are seldom made, due to the difficulty of obtaining mass flow rates to support such measurements. Yet such measurements would provide important information that can be used to help evaluate the total thermal heat transport associated with the Yellowstone Caldera. Typically, discharge from thermal springs migrates through the shallow subsurface, making accurate measurement problematic. Here we present direct measurements of mass and thermal discharge from hot springs in the Lower Geyser Basin of Yellowstone National Park, USA. We added small amounts of nearly pure D2O to four springs in the Morning Mist Springs area that ranged in temperature from 74 to 95 °C and analyzed time-series δD samples to determine the volumes and discharge rates of the test springs. D2O was chosen to limit the ecological and/or visual impacts of other common tracers, such as NaCl or fluorescein dyes. We calculated spring volumes to range between 560 and 27,400 L and estimated mass and heat discharge as 0.08–1.25 L/s and 0.0189–0.312 MW, respectively. The volumes calculated by deuterium doping were larger in every case than those estimated by field inspection, suggesting that the volume participating in shallow fluid circulation is generally larger than is apparent from the surface. The heat flow data, when paired with conductive heat loss estimates in the vicinity of the springs, suggest that current estimates of thermal discharge at Yellowstone may underestimate heat loss from the caldera and offer insights on the rate of magma supplied by the mantle. Thermal flux estimates suggest that a minimum of 3.2–6.3 km3 × 10–2 of basalt magma enters the base of the crust annually.


Hydrothermal heat and mass flux are important components of the Yellowstone Caldera magmatism (Lowenstern and Hurwitz, 2008; Vaughan et al., 2012; Hurwitz et al., 2012). Most of the heat leaving the surface is ultimately derived from the heat entering the base of the crust via both conductive flux across the mantle–crust boundary and as advective heat flux carried across the contact by molten basalt and perhaps by volatiles such as CO2. A better definition of the surficial heat flux may yield a more complete understanding of the material balance that characterizes the volcanic system and, by analogy, older rhyolitic volcanic centers along the Snake River Plain/Yellowstone igneous track (Pierce and Morgan, 1992). The fluxes also play a significant role in determining the caldera’s eruptive potential (Lowenstern and Hurwitz, 2008), and inflation/deflation of the caldera may correlate with changes in heat flux or subsurface hydrothermal fluid migration (Fournier, 1989; Dzurisin et al., 1990).

The primary method for estimating surface heat flow at Yellowstone has been to monitor chloride flux in the streams draining the thermal areas (Fournier et al., 1976; Fournier, 1989; Friedman and Norton, 2007; Hurwitz et al., 2007). This technique assumes that all the discharged fluids were derived from a deeper reservoir with a uniform chloride concentration and temperature. The stream-transported chloride is used to back calculate the total advective heat flow through the surface. Chloride flux is a useful method for estimating energy flow from the Yellowstone Caldera, but it is subject to some assumptions and errors that introduce uncertainty into the results. For example, chloride may be added to the rising fluids via water–rock exchange with rhyolites in the subsurface (Lowenstern and Hurwitz, 2008); furthermore, the proportion of discharge that is carried advectively through the walls of the flow paths feeding hot springs, and the fate of that fluid, are not well known. Some could discharge into streams, and some could recharge the hydrothermal system. It is also not clear the extent to which the chloride flux method accounts for the fraction of conductive heat flow lost through the land surface (Ellis and Wilson, 1955; Vaughan et al., 2012; Hurwitz et al., 2012) in the overall evaluation of the heat budget. Vaughan et al. (2014) estimated that radiant geothermal heat flux was ∼30–45 percent of the total advective heat flux from the hot springs.

Estimating the total thermal and mass flow from the Yellowstone Caldera should rely on integrating flows measured at individual springs. However, there are thousands of identified thermal features at Yellowstone, so to directly measure all of them is not reasonable. Measuring the physical discharge from a spring can be difficult because some or even all of the discharge may occur through the conduit walls below the surface (Benseman, 1959), and the proportion of discharge via steam or other volatile loss may be large (especially in acid sulfate springs). Benseman (1959) pioneered the estimation of hydrothermal spring discharge by adding a conservative tracer to the spring and measuring its decay in concentration over time as a way to overcome the difficulties of directly measuring spring discharge. However, the application of this tracer-decay method to the Yellowstone springs has limitations. Many useful tracers such as NaCl or radioisotopes may be disruptive to ecosystems, particularly where the springs are reservoirs for biodiversity, or they may have unacceptable visual impacts such as with dyes. To overcome these problems, in the present work small volumes of pure D2O were substituted for the more-disruptive tracers. Aliquots of D2O were added to several hot springs in the Mary Mountain Springs area of Lower Geyser Basin, Yellowstone Caldera (Figs. 1 and 2), and then time-series samples were collected from the target springs. The sample δD values were later measured in the laboratory to provide time-decay curves, from which fluid recharge rates were estimated for the springs. The spring discharge estimates reported here were undertaken in conjunction with a program of high-resolution shallow ground temperature measurement (Lubenow et al., 2016) as part of wider study to examine conductive/advective heat partitioning in the near-fields of individual hydrothermal springs.


More than 10,000 hydrothermal features, including geysers, mud pots, fumaroles, and thermal springs, are scattered widely within and around the Yellowstone Caldera (Hurwitz and Lowenstern, 2014) and are distributed over an area larger than 6000 km2 (Fig. 1) (Allen and Day, 1935; White et al., 1971; Fournier, 1989; Christiansen, 2001). This hydrothermal activity is an integral part of the evolution of the Yellowstone Caldera and the associated Yellowstone Plateau Volcanic Field (YPVF). The hydrothermal system has been and continues to be more or less contemporaneous with the caldera magmatism (Larson et al., 2009), and the caldera’s permeable geologic features such as faults or contacts between volcanic units control the distribution of the thermal springs.

The high-silica rhyolite Lava Creek Tuff erupted from the Yellowstone Caldera at ca. 630 ka (Lanphere et al., 2002; Matthews et al., 2015). Since then, ∼25 high-silica rhyolites are known to have erupted within the caldera, mostly as domes and flows, with ages that range from 479 to 70 ka (Christiansen, 2001). Basalts and extracaldera rhyolites that are coeval with the intracaldera rhyolites erupted outside the ring fractures (Christiansen, 2001). No post-collapse basalts are found within the caldera. The magmatism is still active, and seismic tomography has revealed a magma reservoir beneath the caldera (Farrell et al., 2014). The imaged magma body is 90 km long and 5–17 km deep, and it extends 15 km northeast of the caldera. Radiogenic and stable isotope analyses show that all of the rhyolite magmas contain significant crustal components but that the basalts are mantle-derived and are most likely the primary energy source for melting and/or assimilating the crust to produce the rhyolite magmas (e.g., Hildreth et al., 1991; Bindeman and Valley, 2001; Pritchard and Larson, 2012; Pritchard et al., 2013). The Yellowstone basalt and rhyolite magmas are the youngest part of the Miocene and younger Snake River Plain/Yellowstone volcanic progression (Pierce and Morgan, 1992; Leeman et al., 2008).

Most of the thermal springs in the region are concentrated along structural features or along contacts between flows; for example, they are coincident with the presumed locations of the caldera ring faults, along basin–range structures extending north of the caldera, or around the two resurgent domes within the caldera (Fig. 1) (White et al., 1971; Christiansen, 2001; Larson et al., 2009; Hurwitz and Lowenstern, 2014). Measured alteration ages range from ca. 400 ka to younger than 3 ka (Sturchio et al., 1994; Larson et al., 2009). The known ages of hydrothermal activity do not correlate exactly with the timing of post-collapse intracaldera or extracaldera magmatism (Larson et al., 2009). However, this could be an artifact of the dearth of hydrothermal ages older than ca. 100 ka. Younger eruptions may cover older alteration, and extensive glaciation and hydrothermal overprinting may have removed or obscured much of the older evidence. It is not known if the hydrothermal activity has been sporadic since caldera collapse, or, if continuous, if it has waxed and waned in intensity.

White et al. (1971) reported that scientific drilling at Mud Volcano in Yellowstone Caldera had intersected a subsurface steam reservoir. They developed a vapor-dominated model for Yellowstone that incorporates data from vapor-dominated systems at Lardarello, Italy, and at The Geysers, California. This model requires a high heat supply, which is provided at Yellowstone by the 5- to 17-km-deep magma (Farrell et al., 2014). In this model, a low-permeability cap rock is an important component of the vapor-dominated system because it restricts the amount of recharge and thereby helps maintain the system (Ingebritsen and Sorey, 1988). Within the vapor-dominated reservoir, convective steam upflow is concentrated in the more permeable zones such as faults, and downflow of condensate occurs in smaller pores and fractures. Little of the steam escapes to the surface, as most condenses in the shallow part of the system, where it encounters cooler shallow ground water derived from local precipitation.

Both acid sulfate (vapor-dominated) and alkali chloride circumneutral (hot-water-dominated) fluids discharge from Yellowstone’s active thermal springs, commonly in close spatial association (Allen and Day, 1935; White et al., 1971, 1988; Fournier, 1989; Lowenstern and Hurwitz, 2008). Mineral assemblages produced during subsurface hydrothermal water/rock interaction by each fluid type are exposed and in contact with each other in the walls of the Grand Canyon of the Yellowstone River (Larson et al., 2009). Here, it appears that the acid sulfate assemblage overlies the circumneutral assemblage. A genetic link between the two fluid chemistries seems obvious, but the nature of that relationship remains enigmatic. Fournier (1989) proposed that both fluids are derived from a single deep parent fluid at a temperature of ∼340 °C and containing ∼400 ppm chloride ion and abundant dissolved CO2 and H2S.


The Lower Geyser Basin (LGB) is a large hydrothermally active area that lies north and west of the Mallard Lake resurgent dome (Fig. 1). Four springs were chosen for study in the northwestern part of the LGB within the drainage of Nez Perce Creek (Fig. 2). The springs were selected because they fit several criteria for the experiments: (1) All four springs have diameters less than a few meters, and therefore their volumes allow adding manageable amounts of D2O with ease of promoting mixing to raise their δD values. (2) The springs are individually isolated and are at least tens of meters away from other larger springs so that overlapping thermal effects could be avoided as much as possible. (3) All four springs have average temperatures higher than 70 °C and appear to have vigorous recharge, although one of them (LCBNN159; Rodman and Guiles, 2008, “Bee”) has no visible surface discharge. (4) All four springs are easily accessible via a short hike from the Grand Loop Road (U.S. Highway 89/191) on the Mary Mountain/Nez Perce Creek Trail and are shielded by forests and topography from view by park visitors, thus allowing both ease of access and minimal visual impact. (5) All four springs lack visible evidence of biotic communities on their walls (i.e., no microbial mats). Finally, (6) all the research tasks could be performed in a safe manner without negatively impacting the springs or their environment. Table 1 lists the four springs with their locations and other characteristics. Informal names were assigned in the field (Aspen, Bee, Carolyn, and Edgar) and are used in this study. Formal spring names and Yellowstone National Park Research Coordination Network designations are found in Table 1 (Rodman and Guiles, 2008).

The study area (Fig. 2) is in the Morning Mist Springs area, where a sparsely vegetated meadow is underlain by hot-spring sinter deposits and pockmarked by ∼160 vigorous to inconspicuous thermal features (Muffler et al., 1982). It lies between Nez Perce Creek to the north, the Elephant Back rhyolite flow to the south and southeast, and the glacial kame deposits (Porcupine Hills), which are locally cemented by hydrothermal minerals to the west (Christiansen, 2001). Nez Perce Creek flows along the front of the Nez Perce Creek rhyolite flow (148 ka), which overlies the Elephant Back flow (153 ka) beneath the study area (Muffler et al., 1982). Two U.S. Geological Survey research drill holes, Y-4 and Y-13, were sited in the Morning Mist Springs area (Fig. 2) (White et al., 1975). Y-13 was drilled within a few meters of the Porcupine Hills Geyser (Aspen) and was placed to intersect the center of a hot, upflow zone, whereas the Y-4 site was chosen because it was some distance away from the active thermal features. Similarity between the fluids from both wells led White et al. (1975) to suggest that a single deep hot fluid reservoir was their source.

Field measurements and experiments were carried out between June 22 and June 28, 2014. The springs’ surface areas and depths were measured in the field, and their volumes were calculated assuming ideal cone or cylinder shapes (Table 1). Complications in estimating volume were encountered at Aspen, because it was deeper (>6 m) than the available sounding measure, and at Edgar, because it was too wide to allow access to measuring its total depth at its center. Springs were located with handheld GPS units. Temperatures were measured at a variety of depths in each spring using a digital thermometer and K-type thermocouple. Higher spring temperatures were desired because they suppress biotic communities, which could interfere with mixing and could perhaps partition D2O inequitably and impact the results. Relatively constant high temperatures also indicate that a spring is well mixed. Ambient air temperatures were monitored at 45-second intervals for the duration of the study.

Once spring volumes were estimated, D2O volumes appropriate to produce an estimated change in the springs’ δD values of ∼150‰ were measured using a graduated cylinder (Table 1). Laboratory comparisons between the graduated cylinder readings and pipetted volumes shows the precision of the spike volume to be about ± 0.25 ml. Once analyses were conducted in the laboratory, it became evident that the springs volumes as estimated in the field were too small and that the shifts in the springs δD values were in the range 50–100‰. The heavy water was transferred to a lab-grade Nalgene bottle. The D2O was carefully poured into the spring at time t = 0. The bottle was then washed quickly with spring water and emptied back into the spring to ensure complete D2O transfer. A collapsible oar was used immediately after t = 0 to carefully stir the spring to facilitate mixing. Care was taken to avoid disturbing any sediment with the oar. Three stopwatches were started synchronously at t = 0 so that the sample collection times after t = 0 were accurately recorded. At least one experiment was conducted on each of the four springs; however, Caroline was tested twice (Carolyn 1 and Carolyn 2).

A 1.5 m × 2 cm polyethylene pipe with a single-hole rubber stopper in one end was used to collect the samples. The end with the stopper was lowered into the spring, and the sample was extracted by covering the other end of the pipe in much the same way a drinking straw can be used to pipette a liquid. The samples were collected from the same central point within each spring, usually at one half the depth of the spring. The sampling tube was marked on its side for each spring so that the sample depth was consistent. Samples were stored in centrifuge tubes sealed with paraffin and were placed out of light for transport back to the laboratory. No reagents were added to the samples in the field. Several background samples were collected from each spring prior to introducing the heavy water. A sample was collected closely after t = 0 and then at shorter (several-minute) intervals for about an hour. After the initial hour, the intervals between sampling were increased, such that between 25 and 30 samples were collected over a period of 6 h to nearly a day, depending on the volume of the spring.

Analyses were conducted in the stable-isotope laboratory in the GeoAnalytical Laboratory at Washington State University using a ThermoFinnigan Delta V+ gas-source isotope ratio mass spectrometer with a ConFlo IV interface. δD values were measured using pyrolysis with a high-temperature conversion elemental analyzer (TC/EA). Using an autosampler, ∼500-nL volumes were drawn into a 1.2-µL syringe and were injected directly into the TC/EA unit at 1400 °C. The resultant gas was carried by a continuous flow of He to the mass spectrometer. Samples were analyzed sequentially for 5 to 12 times to overcome memory effects. The δD values of all the samples were measured.

δ18O values were measured by H2O equilibration with CO2 using a ThermoFinnigan GasBench II. After 24 h of equilibration at 32 °C, the autosampler was used to transfer the CO2 to the mass spectrometer. Values were corrected for fractionation between CO2 and H2O using the factors of Bottinga (1968). The δ18O values of all the background and selected tracer-experiment samples were measured. Replicate analyses of standards show that the δD and δ18O measurements are precise to ± 1.0‰ and ± 0.1‰, respectively.

During the course of the study, three samples of Yellowstone meteoric water were collected. Two rain samples were collected during the tracer experiments, one on June 24 and one on June 26, 2014. The June 24 sample contained some pea-sized hail, whereas the June 26 sample was pure liquid. The third sample was collected from a snow berm along the Grand Loop Road near Mt. Washburn at latitude 44.7757°N and longitude 110.4596°W. The δD and δ18O values of these three meteoric water samples were measured together with the hot springs samples. They yielded O ratios of –14.72‰, –14.04‰, and –16.83‰, and H ratios of –108.67‰, –104.94‰, and –123.62‰, all relative to international water standard VSMOW (Vienna standard mean ocean water), respectively.

All the measured δD and δ18O values were normalized relative to the international water standards VSMOW2, GISP, and VSLAP2, and relative to the internal standards DSW (desalinated seawater) and WAWA2 (Washington water, which is deionized local tapwater), all of which were analyzed along with the samples. All the unknown sample analyses were bracketed by the standard analyses, and normalization was completed using the method of Nelson (2000). A normalization calibration curve was constructed by plotting the measured values versus the accepted values of the standards. The curve was applied to the measured sample values. The isotope ratios and other data for the tracer-decay experiments are shown in Tables 26 and plotted in Figures 47. All the isotope ratios are reported in the familiar delta notation relative to VSMOW.


Background and Meteoric Water Isotope Ratios

Background samples were collected because their analyses define the baseline for δD ratios in the hot springs to which the tracer-decay experiments have to be calibrated. Their 18O/16O ratios complement the D/H ratios and place the test spring waters within the context of the LGB springs and the entire Yellowstone hydrothermal system. The background stable-isotope ratios for the Morning Mist Springs samples plot near or within the field of other LGB hot springs fluids (Fig. 3) (Ball et al., 2006, 2010; Holloway et al., 2011).

The three precipitation samples analyzed in this study plot near or along the meteoric water line (Fig. 3) and within the range of Yellowstone meteoric water ratios (Thordsen et al., 1992; Rye and Truesdell, 2007). The Local Meteoric Water Line (LMWL) at Yellowstone, δD = 7.97δ18O + 9.67 (Thordsen et al., 1992), lies very close to the meteoric water line, δD = 8δ18O + 10, as originally defined by Craig (1961) (Rye and Truesdell, 2007). Kharaka et al. (2002) found that precipitation samples from a high-elevation LMWL at Yellowstone lie along the line δD = 8.2δ18O + 14.7. The array of all the LGB hot springs fluids project back to the meteoric water line near the hydrothermal-recharge parent fluid of Truesdell et al. (1977) at δD = –149‰ and δ18O = –19.9‰. The parent hydrothermal fluid has lower isotope ratios than the modern meteoric water samples from this study and Thordsen et al. (1992).

Fluids from other hot springs basins at Yellowstone also lie along the LGB trend in Figure 3, or broadly along its linear extension (Truesdell et al., 1977; Ball et al., 2006, 2010; Rye and Truesdell, 2007; Holloway et al., 2011). The hot springs fluids throughout the region are derived mostly from a deeper isotopically homogeneous fluid with little addition of shallow modern groundwater. The fluid could be an older meteoric water that recharged during a cooler climate in Pleistocene time (Rye and Truesdell, 2007), as has been demonstrated for other aquifers in the Pacific Northwest where deeper older waters have lower δ18O ratios than shallower younger water (Larson et al., 2005). Alternatively, the aquifer was recharged from higher peaks around the caldera (Kharaka et al., 2002; Rye and Truesdell, 2007), where higher elevation precipitation tends to have lower isotope ratios than at lower elevations (Alley and Cuffey, 2001). The parent hydrothermal fluid must have experienced deeper boiling prior to rising in the vapor-dominated hydrothermal system. The isotope data suggest that the volume of the rising fluid most likely dwarfs any significant mass contribution from the local groundwater.

Hot Springs Volumes

The hot springs volumes were first estimated in the field (“Field Volumes” in Table 1) to determine how much tracer to add to shift the δD values to provide a sufficiently large signal to accurately track its decay. These calculations were used to estimate the physical volume of the spring discharge pools, which, for convenience, are referred to as “cisterns.” Although these calculated volumes were useful for estimating the amount of tracer to add to sufficiently perturb a given spring, the cistern volume may not represent the volume of the spring water that mixes with the D2O. There is an interface between the inflow channel and the cistern volume that accommodates mixing, and that interface most likely does not lie where the channel and cistern exactly intersect; instead, incoming fluids probably interact with fluids resident in the cistern in a complex way across a poorly defined mixing volume. For example, the spring fluids cool at the surface and the cooler fluid will tend to flow downward along the sides of the cistern, reentering and mixing with the inflowing fluids in a buoyant convection roll. This was observed in several of the springs, most notably in Carolyn. Therefore, the actual volume of the spring where mixing between spring water and the added D2O takes place is designated as an effective mixing volume (EMV). This volume may be similar to but not identical with the physical volume defined by the geometrically idealized shape of the cistern.

Although the idealized-geometry estimate of the cistern volume is useful for estimating D2O additions, it is not sufficiently accurate for the estimation of recharge/discharge rates or advective heat fluxes. A more precise estimate of the EMV for a given spring cistern can be obtained using mass balance. Because the change in the δD value of the hot spring between background and the “spiked” value at time t = 0 is a function of the amount of pure D2O added and the volume into which it is mixed, the EMV can be calculated from the change in D2O concentration (before and after spiking) and the number of moles of deuterium added: 
where ΔC is the change in the concentration of deuterium [mol/L] when A moles of D2O are added to a spring with a volume V [L]. If the volume change from the addition of A moles of D2O, which is four to five orders of magnitude less than the cistern volume, is assumed to be negligible, the EMV (as V) is given by: 

The number of moles, A, of pure D2O added to the spring can be calculated from the volume of D2O measured in the field (Table 1) and the molarity of pure D2O as 55.27 mol/L. It should be noted that the equation for calculating a δD value includes a term for the ratio D/H; however, pure D2O has no H, and its δ value is undefined. As a result, D2O molarities were calculated from the δ equation using the measured background δD values and a D/H ratio for VSMOW of 1.558 × 10–4. The precision of A, the measured volume of added heavy water, is better than 1%, and the H isotope ratio is precise to within 1‰. This suggests that the precisions of the calculated EMVs are near 1%.

Recharge/Discharge Rate Calculations

If the volume of a reservoir remains constant, then by mass balance the recharge and discharge rates can be considered equal as long as all input and outflow terms are accounted for (i.e., discharge through a well-defined channel, losses through the walls of the spring conduit, etc.). The rate at which D2O is removed from the reservoir is proportional to the amount of D2O remaining in the spring; therefore, the return to background can be modeled as an exponential decay process as a completely mixed flow reactor. Assuming constant cistern volume for the time period of monitoring, the model is given by: 
where C(t) is the concentration of D2O in the spring [mol/L] at time t [s], C0 and Ci are the initial (spiked) and background concentrations in the spring, respectively, and tc is the characteristic time of the spring [s]. The characteristic time is given by: 
where V is the EMV [L] and Q is the recharge/discharge rate [L/s], which are assumed to be equal. The characteristic time is the time at which the concentration of D2O in the spring will have dropped to 1/e (∼37%) of the initial concentration. This characteristic quantity is often referred to as the “residence time” of a reservoir. The residence, or characteristic, time can be determined for a given spring by adjusting the value of tc until the exponential decay plot matches the spring data for the return of the spring to background values of D2O. Because the change between background and the initial concentration, (C0Ci) = ΔC in Equations (2) and (3), is known from the data, the EMV (as V) can be found from Equation (2), and the modeled characteristic time then can be used to calculate the value of Q using Equation (4). The data from each of the five experiments are shown in Tables 2 through 6 and Figures 4 through 7, and the spring characteristics, including the estimated recharge/discharge rates (Q), are provided in Table 7. The recharge/discharge rates should be precise within a percent, as the δD values are ± 1‰ and the EMVs are ± 1%.

Recharge/Discharge Rates

Recharge/discharge rates for the hot springs in Table 7 range over about two orders of magnitude, from 0.08 L/s for Bee (LCBNN159) to 1.25 L/s for Edgar (LCBNN161). Similarly, the characteristic, or residence, times of the four springs range from 1.9 h for Bee to 15.4 h for the second of two tests performed on Carolyn.

One point worth noting is the difference between the springs (cisterns) volumes estimated by approximating the springs as idealized geometric shapes and the calculated EMV volumes. In all cases, the EMVs exceed the volumes estimated from the idealized geometric shapes. One particular example of this is the volume of Aspen (Porcupine Hills Geyser), which is cylindrical in shape to the extent visible from the surface. However, obtaining accurate measurements of the depth of the spring cistern proved to be problematic due to the large apparent depth, and it is possible the true depth could be significantly greater than the six or so meters that were measured in the field. Note that Aspen’s EMV (Table 7) is significantly larger than its field volume (Table 1).

The disparity between the estimated and calculated volumes for spring Carolyn was also much larger than anticipated. The field estimate for the cistern volume, made on the basis of surface observations, was 1532 L (Table 1); however, data from the first Carolyn tracer test (Carolyn 1) place the EMV at 9173 L. Although the volumes calculated from mass balance were expected to be larger than the estimates made from surface observations, the magnitude of this discrepancy for Carolyn was surprising and suggests that the D2O tracer was recirculating and likely mixing with waters at greater depths and with a much larger reservoir than the volume of the spring cistern visible from the surface. As discussed above, buoyant convection currents (upwelling hot fluids and downwelling cooler flows) were seen to descend along the outer edges of the cistern, apparently reentering the conduit by which rising hydrothermal fluids were discharging into the cistern from depth. It is probable that there is a large volume concealed below and that it is connected to the Carolyn cistern by this inflow conduit.

The unexpectedly large EMV for Carolyn somewhat complicated tracer-test analyses for that spring. Two tracer-decay experiments were performed on Carolyn, with the objective of checking test reproducibility. The two tests, designated Carolyn 1 and Carolyn 2, were carried out on June 22, 2014, and June 28, 2014, respectively. The two tests were conducted ∼110 h apart with the expectation that enough time had elapsed for the δD value to return to background. However, subsequent analyses on the two “background” samples collected prior to the Carolyn 2 tests showed that the second experiment began with a mean value of –132.6‰ for δD, somewhat higher than the –139.7‰ value determined for the Carolyn 1 test. This shift in the apparent background implies that ∼9.9 mL of the original 70 mL “spike” of D2O was still present in the EMV at the start of the second test. To compensate for the elevated initial condition, the spike volume for the Carolyn 2 test was calculated for 79.9 mL added D2O, instead of the 70 mL added at t = 0. With this correction, the characteristics of the exponential decay δD versus time plots for the two Carolyn experiments are similar. However, the fact that the background δD value had not been reestablished more than 110 h after the first of the two experiments provides additional evidence to support the conclusion from the mass balance/EMV calculation that the mixing volume is significantly greater than the cistern volume visible from the surface, and this may suggest even more complicated plumbing issues deeper in the system conduit.

Heat Transfer

Fluids discharging from the hot springs also carry heat energy by advection. Once recharge/discharge rates for the individual springs are known, the quantity of heat transferred by advection can be calculated by: 
where qadv,out is the advection of heat through hydrothermal water, Q is the spring discharge calculated using equation (4), ρ is the density of discharging water (999.97 kg/m3), c is the specific heat of discharging water (4186 J/kg K), and Ts is the temperature of the discharging water as it leaves the spring (assumed to be the well-mixed spring temperature) (Table 7). As a rule, the rate of heat transfer is not an absolute value but is relative to some reference temperature, Tr in Equation (5), taken here to be the average ambient air temperature of 11 °C (Ingebritsen et al., 1989). Table 7 shows the qadv,out for each experiment calculated using equation (5).

Given that the heat transfer for any individual spring is dependent on its discharge temperatures and on the discharge rate calculated using Equation (5), a comparison of the advective heat flux can be made between hydrothermal features if the same reference temperature is used for all features. Unfortunately, there are few individual hydrothermal features at Yellowstone that have modern discharge rate estimates that can be used for comparison to the results presented here, although temperature data are abundant in the literature. Excelsior Geyser in Midway Geyser Basin is one feature that does have a surface discharge rate estimate. According to the National Park Service, the geyser has a discharge rate of 4000 gal/min (255.5 L/s) during noneruptive periods and a water temperature of ∼93 °C. Heat output is estimated to be 8.5 × 107 W using Equation (5) and the same values for the other parameters used here. The total heat output of the hot springs included in this study is only 5 × 105 W. This disparity is expected given the orders-of-magnitude differences in discharge rates between Excelsior and the four smaller springs discussed here.

The chloride inventory method estimates heat flow carried by liquid water discharging from the thermal springs to range from 4.5 to 6.0 GW (Fournier, 1989; Friedman and Norton, 2007). The chloride flux data yield a total hot springs liquid discharge rate of between 3090 and 3960 L/s for the Yellowstone hydrothermal system (Fournier, 1989; Hurwitz and Lowenstern, 2014; Eppelbaum et al., 2014). It is interesting to note that peaks in chloride flux occur annually during the spring months coincident with snowmelt runoff (Norton and Friedman, 1991). The spikes have been hypothesized to result from increased thermal discharge due to a snowmelt-driven increased hydrostatic head (Friedman and Norton, 1990), release of excess chloride stored in lakes and the shallow subsurface (Ingebritsen et al., 2001), or higher discharge from groundwater due to a higher hydrostatic head during periods of heavy rainfall (Friedman and Norton, 2007). Seasonal variations are therefore important for understanding Yellowstone hot springs discharge and advective heat transfer. Although the results presented here do not account for these seasonal effects, the method could be adapted to include those influences.

Thermal Flux from the Yellowstone Plateau Volcanic Field

The enhanced heat flow across the Yellowstone surface is ultimately derived from heat transferred across the mantle–crust boundary, both from conduction and by magmatic advection and other fluid flow. If accurate estimates of surface heat discharge can be obtained, they can be used to place limits on the heat input into the base of the crust. Much of the work on determining an overall heat budget for the YPVF has focused on obtaining and analyzing the flux of chloride in surface waters, under the assumption that all the chloride leaving the YPVF originated from a single parent fluid having a chloride concentration of around 400 mg/L and a temperature of 340 °C (Fournier, 1989; Friedman and Norton, 2007). The exact chloride concentration, homogeneity, and temperature of the parent fluid are debatable, however, and Hurwitz and Lowenstern (2014) used various assumed initial compositions and temperatures to estimate total heat flux values between 4.0 and 8.0 GW. Estimates based on chloride flux are closely tied to the advective flow of hydrothermal fluids from hot springs and do not account for conductive heat transport across the surface boundary. Chloride flux monitoring, therefore, provides a minimum estimate on heat transfer across the mantle–crust boundary.

Estimations of heat flux in the YPVF play a critical role in understanding deep processes and the material balance of the Yellowstone volcanic system. The amount of heat supplied by solidifying basalt magma that has moved into the crust across the basal crustal boundary depends on the composition of the magma, how much it cools, and where it solidifies in the crust. This energy/material balance is demonstrated here, making a simple calculation using the software rhyolite-MELTS (Gualda et al., 2012; Ghiorso and Gualda, 2015) to examine the enthalpy produced by typical basalt magmas as they change from liquidus to solidus states. An ordinary parental basalt (∼8.5 wt% MgO) with 1 wt% H2O at 11 kbar (equivalent to 40-km-deep crust at mean density of 2800 kg/m3) is fully liquid at 1335 °C with an enthalpy of –11.7 kJ/g. MELTS predicts the basalt to be almost entirely crystalline (4 wt% liquid) at 700 °C with a total (melt and crystals) enthalpy of –12.9 kJ/g, with a net difference of 1.24 kJ/g released over the 635 °C drop in temperature. The same calculation for the magma at 7 kbar (∼25 km depth), where the liquidus is near 1270 °C, shows a net enthalpy release of 1.11 kJ/g. In both cases, this yields ∼1.95 J/gK.

A picritic parent may be more appropriate for Yellowstone. Using the picritic glass from Kilauea, Hawaii, at ∼16 wt% MgO and 0.34 wt% H2O (Clague et al., 1995; Norman and Garcia, 1999), MELTS puts the liquidus at 1517 °C at 11 kbar. The enthalpy change to 700 °C is 1.53 kJ/g, which yields 1.88 J/gK over the 817 °C drop in temperature. For all cases examined here, the enthalpy release is ∼1.9 J/gK, and this value is used in the thermal balance model.

For this model, the temperature drop in the magma is defined by the change in temperature from a liquidus at 1400 °C to a solidus at 700 °C. These temperatures are based loosely on the MELTS calculations. Thus, over the 700 °C drop in temperature at 1.9 J/gK, the enthalpy released is 1330 J/g. The heat flux estimated by the chloride flux model (Hurwitz and Lowenstern, 2014) yields a total advective heat transfer of 4.0–8.0 × 109 J/sec. Dividing this advective flux by the magma enthalpy release yields 3.01–6.02 × 106 g/s of cooled magma. Thus, a simple thermal balance suggests significant and ongoing basalt magma influx into the base of the crust under Yellowstone. Although this is a simplistic calculation that ignores other crustal heat sources or sinks, some of which are most likely not trivial, it does serve as a first-pass assessment. Assuming a basalt magma density of 3 g/cm3, this translates to 1.00–2.01 × 106 cm3/sec, or 3.2–6.3 × 10–2 km3/yr.

In comparison, the long-term average eruption rates at Kilauea are ∼0.1 km3/yr (Swanson, 1972; Poland et al., 2012). The chloride flux-based thermal model, therefore, yields basalt input into the Yellowstone crust at rates one order of magnitude less than Kilauea eruption rates. Lowenstern and Hurwitz (2008) used CO2 flux at Yellowstone to suggest a basalt influx into the lower crust of ∼0.3 km3/year. This estimate is an order of magnitude greater than that presented here, and three times the Kilauea eruption rate.

Lubenow et al. (2016) estimated near-field conductive heat losses around Bee (LCBNN159) on the basis of 20-cm-deep temperature measurements in a grid at a 3 × 3 m spacing around the spring vent. Their results suggest that in some cases, conduction may account for as much as 50 percent of the total heat transfer to the atmosphere across the land surface boundary, where the remaining flux of heat, calculated here, is included in the advective spring flow. If a 50:50 (conduction:advection) distribution of heat transfer is general to the Yellowstone hot springs, then the heat fluxes estimated from stream monitoring of chloride flux may account for only half the total Yellowstone heat discharge, and the estimate of the basalt transport across the mantle–crust boundary would be at least twice the 3.2–6.3 × 10–2 km3/yr calculated from the heat flux estimates based solely on the chloride balance method.


The experiments described here have shown that D2O can be used as a reliable tracer in alkaline-chloride hot springs. The tracer-decay experiments demonstrate that the discharge of hydrothermal fluids from the springs can be modeled as a completely mixed flow reactor. Together with the ability to measure tracer concentrations precisely, precise volume, discharge rate, and advective heat transfer estimates can be calculated.

The use of labeled water in the tracer-decay experiments at Yellowstone poses no threat to the environment within and around the hot pools. The addition of the D2O to the hot springs, in fact, did not increase their δD values above that of seawater for most of the duration of the experiments, while providing a strong signal, the decay of which was easily measured and modeled. The characteristics of D2O make it a nearly ideal tracer for use in hydrothermal springs: It does not disturb or harm any of the organisms in their natural environment, it is ephemeral and dilutes quickly, it is not visible, and it has no adverse effects on the geochemistry of the springs.

Each of the experiments returned results that could be fit to the completely mixed flow reactor model. Mixing took place over longer time periods than expected (up to ∼5400 seconds), but initial concentrations could be back-calculated by using only data collected after mixing had been completed and the linear nature of discharge on a log plot had been established. Volumes (EMVs) for each of the springs were calculated from the initial (spiked) concentrations; these volumes ranged from 563 L (Bee) to 27,399 L (Edgar). Simple field measurements of spring volumes do not accurately provide a measure of the true volumes (the EMVs) into which the tracer mixes. Discharge rates for the hot springs based on the EMVs ranged from 0.08 to 1.25 L/s (Bee and Edgar, respectively).

Data obtained from the two experiments at Carolyn yielded almost identical decay curves. Because residual D2O tracer from the first experiment was still present at the beginning of the second, values were slightly higher in the Carolyn 2 experiment. When this excess tracer was accounted for in the calculations for Carolyn 2, it yielded results almost identical to those from Carolyn 1. The slow, long-term decay in Carolyn suggests that there is a concealed reservoir beneath its cistern that communicates with it. It is possible that D2O tracer-decay experiments can provide information about deeper plumbing in a hot spring, and future work is planned to address this potential application of the method.

The chloride flux model appears to provide a minimum estimate of the total thermal discharge from the YPVF. Uncertainties in the homogeneity of the hydrothermal discharge across the entire volcanic field contribute to underestimating discharge. However, the most significant factor in the underestimate is most likely that conductive heat transport across the surface was not included. The chloride flux estimates of 4.0–8.0 GW for Yellowstone can be balanced by an input of 3.2–6.3 × 10–2 km3/yr of basalt magma advectively transporting heat into the base of the crust. There are many uncertainties in this calculation, and this basalt flux should be considered a minimum estimate.


Many people assisted this research at various stages. We sincerely thank the rangers in the research permit office at Yellowstone National Park for their interest, support, and encouragement. Charles Knaack in the GeoAnalytical Laboratory at Washington State University provided invaluable technical assistance for the analyses. Many have also contributed to the fieldwork and are gratefully acknowledged, including Keegan Schmidt, Brady Lubenow, Alayne Donnelly, Gilbert Ching, Jay Meyers, Ben Jones, Jennifer Light, Alex Moody, and Erika Rader. Steve Nelson and Jake Lowenstern provided useful and stimulating reviews of the manuscript, together with useful comments from Eric Christiansen. Tom Sisson provided an introduction to the rhyolite-MELTS calculations and some very useful discussions, for which he is gratefully thanked. The research was supported by funding from the National Science Foundation (NSF1250435, PI: Larson; NSF1250381, PI: Fairley). The work described was performed under permit from the National Park Service (permit YELL-2014-SCI-6034).

Science Editor: Shanaka de Silva
Associate Editor: Eric H. Christiansen
Gold Open Access: This paper is published under the terms of the CC-BY license.