References

Abstract. A unique conceptual model describing the conductive heating of rocks in the thick unsaturated zone of Yucca Mountain, Nevada by a silicic pluton emplaced several kilometers away is accepted by the US Department of Energy (DOE) as an explanation of the elevated depositional temperatures measured in fluid inclusions in secondary fluorite and calcite. Acceptance of this model allowed the DOE to keep from considering hydrothermal activity in the performance assessment of the proposed high-level nuclear waste disposal facility. The evaluation presented in this paper shows that no computational modeling results have yet produced a satisfactory match with the empirical benchmark data, specifically with age and fluid inclusion data that indicate high temperatures (up to ca. 80 °C) in the unsaturated zone of Yucca Mountain. Auxiliary sub-models complementing the DOE model, as well as observations at a natural analog site, have also been evaluated. Summarily, the model cannot be considered as validated. Due to the lack of validation, the reliance on this model must be discontinued and the appropriateness of decisions which rely on this model must be re-evaluated.


Introduction
Yucca Mountain in southern Nevada was studied by the US Department of Energy (DOE) between 1983 and 2005 as a prospective host for the first facility for the geological disposal of spent nuclear fuel and high-level nuclear waste in the United States. The prolonged period of study at the site culminated in June 2008 when DOE filed its License Application with the US Nuclear Regulatory Commission (DOE, 2008). In 2009, however, the administration of President Obama decided to terminate the Yucca Mountain project, for reasons purportedly unrelated to the safety or technical suitability of the site (GAO, 2011). Accordingly, in 2010 the DOE filed a motion with Nuclear Regulatory Commission to withdraw the license application. Despite the continued uncertainty regarding the legal standing of this decision, over 2009-2011 the repository program was dismantled and effectively shut down.
At Yucca Mountain, nuclear waste was to be emplaced in the thick unsaturated (vadose) hydrogeologic zone. The safety case for the disposal facility rests heavily on the concept that water is and always has been scarce in this zone (i.e., during the last 11.6 myr; DOE, 2008). In the course of the Yucca Mountain site characterization activities, however, mineralogical and geochemical information became available indicating that waters with elevated temperatures (up to 85-90 • C) had moved through the unsaturated zone in the past. A conceptual model proposed by the US Geological Survey (USGS) explained these temperatures as the conductive heating of the unsaturated-zone rocks by a large magma body emplaced in the late Miocene (ca. 11 Ma) under the Timber Mountain caldera, some 7-10 km north of the repository site (Marshall and Whelan, 2000).
The purpose of evaluating the USGS/DOE conceptual model presented in this paper is to establish whether or not the validation of the model can be considered successful; in other words, whether or not the model is viable. The need for this evaluation stems from two facts: (1) no formal evaluations of this model have been published; and (2) despite this, the model has become firmly established in the literature, and was relied upon in the decision making related to high-level nuclear waste disposal in the USA (e.g., Dublyansky, 2007). The evaluation presented in this paper was performed according to a generic protocol for the evaluation of geoscientific models (Grewe et al., 2012).

Paleotemperatures determined from secondary minerals
Volcanic tuffs in the thick unsaturated zone of Yucca Mountain host secondary minerals that are ubiquitous (primarily, calcite and silica). The minerals were initially interpreted as having been deposited from rain waters that had percolated from the surface (Szabo and Kyser, 1990;Paces et al., 2001;Whelan et al., 2001Whelan et al., , 2002. Subsequently, however, it was established on the basis of fluid inclusion studies of minerals collected in the ESF tunnel complex (Exploratory Studies Facility, a 7.8 km long C-shaped tunnel excavated into the host formation of the planned repository) that the minerals were deposited from waters whose temperature reached 85-90 • C Wilson et al., 2003;Whelan et al., 2008). The age of the oldest fluid inclusion temperature (ca. 77 • C) was constrained to 9.4 ± 0.7 Ma by U-Pb dating (Whelan et al., 2008; Table 4). Deposition of the secondary minerals took place at depths of 30 to 300 m from the contemporaneous surface of the mountain. Temperature increases related to a "normal" geothermal gradient would not be expected to exceed ca. 25-28 • C. An important question arose: how to explain the movement of the conspicuously thermal waters through the rock which, according to the accepted geological understanding (e.g., DOE, 2001), had always been a few hundred meters above the water table during the last 11.6 myr?

Thermal history of Yucca Mountain: what is known?
Cooling of tuffs -The rhyolitic tuffs of the Paintbrush Group that comprise most of the unsaturated zone in Yucca Mountain were emplaced during large-scale silicic eruptions 12.7 Ma. The cooling of the 350 m thick sheet of ash-fall tuff from its estimated temperature of emplacement (680-720 • C) to ambient temperatures took about 7000 years (Buesch and Riehle, 2007). The younger Timber Mountain Group tuffs (ca. 11.45 Ma) may have been deposited on top of the already cooled Paintbrush tuff. Subsequently this later layer has been largely eroded away; the maximum estimate of its thickness is 100 m (DOE, 1995). Therefore, any thermal water that circulated through the unsaturated zone of Yucca Mountain after ca. 11.4 Ma cannot be related to the residual heat of the tuffs. The Timber Mountain caldera hydrothermal event -A large silicic magma body was emplaced beneath the Timber Mountain caldera, 7-9 km to the north of Yucca Mountain shortly after the final climactic eruption at 11.45 Ma. This resulted in the development of a large-scale southward-flowing hydrothermal plume. According to Bish and Aronson (1993), the system included an upflow zone in the area of the Claim Canyon caldera, where thermal waters likely discharged at the surface (Fig. 1). Further to the south they affected only deep parts of the rock sequence (under Yucca Mountain -ca. 1000 m and deeper). A pronounced north-south thermal gradient has been noted: ". . . it is apparent that a significant thermal event has occurred in the northern end of Yucca Mountain but has not significantly affected the southern end" (Bish and Aronson, 1993, p. 153). The authors also argued that the upper ca. 1000 m of Yucca Mountain rocks were significantly cooled by infiltration of meteoric water (the "rain curtain" effect). It is thus apparent from Fig. 1 that temperatures of 75-90 • C at the ESF level cannot be attributed to the Timber Mountain caldera hydrothermal event.
Another constraining issue is given by the time frame. The Timber Mountain caldera hydrothermal system was active between ca. 10. 5-11.0 Ma (Bish and Aronson, 1993, p. 159); i.e., the event ended well before the time when the oldest elevated-temperature secondary minerals sampled in the ESF were deposited.

Validating the MICH model through computational modeling
Although "validation" may not be the most appropriate term (e.g., NRC, 2007;Nordstrom, 2012), for the purpose of this evaluation it will be used according to definitions adopted in documents related to the nuclear waste disposal. Validation means "determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model" (AIAA, 1998) or "building confidence that a model adequately represents a real system for a specific purpose" (IAEA, 2003). The DOE administrative procedure, Analyses and Models (AP-3.10Q), stipulates that: ". . . model validation shall consist of comparing analysis results against data acquired from the laboratory, field experiments, natural and man-made analog studies, or other relevant observations" (p. 13). Validation through computational modeling involves the construction of a numerical model and comparison of the computational outcomes with the experimental outcomes. Failure to pass the quantitative comparison test would result in a need to revise the conceptual model; failure to achieve the acceptable agreement through revisions must lead to the abandonment of the conceptual model.
For validation purposes, the outcomes of computational modeling must be compared to a benchmark. For the MICH model, the benchmark consists of the coupled temperatureage data characteristic of a certain distance from the heat source (magma chamber). The empirical data constraining the "thermal history" of the unsaturated zone of Yucca Mountain were obtained by studying secondary minerals in the ESF tunnel complex, as noted. The results are summarized in graphical form in Fig. 2. The graph was first presented in Whelan at al. (2003, Fig. 4) and then reproduced in Whelan et al. (2008, Fig. 8b), where it was used as a benchmark for the computational modeling. A simplified version of this graph was used as a benchmark in Dublyansky and Polyansky (2007).

Evaluation of the benchmark
The age and temperature data used in the benchmark come from three different sources.
Fluid inclusions -Homogenization temperatures obtained from calcite and fluorite coupled with the 235 U / 207 Pb age dates measured in closely associated chalcedony and opal are shown in Fig. 2 by solid symbols. The data points for which only minimum ages are available (triangles and dotted lines) have too-large uncertainties in the time scale and therefore are not useful. In addition, errors of the temperature determinations (typically, 3-7 • C) are not shown on the graph.
Stable isotope calculations -The circles in Fig. 2  were constrained by U-Pb or U-series ages of silica phases, and were interpolated or extrapolated by assuming a constant rate of mineral deposition (Whelan et al., 2008, Sect. 5.2.2). Considering that known variations of the growth rate in the Yucca Mountain samples exceed one order of magnitude (0.5 to 5 mm Ma −1 - Neymark et al., 2002;Paces et al., 2004;up to 10 mm Ma −1 - Wilson et al., 2003, Fig. 8b), the 1 mm error in determining the position of the sample corresponds to up to 0.5 myr. Given the complexity of growth textures observed in the Yucca Mountain samples, age uncertainties up to 1-2 myr can be expected for the stable isotope data points.
The uncertainty in the isotopic temperatures shown in Fig. 2 is, primarily, the uncertainty of the method. Calculation of the temperature from δ 18 O of calcite involves a number of assumptions, namely: (1) that deposition occurred in isotopic equilibrium with regards to the mineralforming solution; (2) that the equation selected, O'Neil et al. (1969), faithfully describes the temperature-dependent fractionation of oxygen isotopes between calcite and water; and (3) that the mineral-forming water has retained a constant δ 18 O = −11 ‰ throughout the last 10 myr. Each of these assumptions has an inherent uncertainty. Assumption 1 -Deposition under isotopic equilibrium may be theoretically compatible with the MICH conceptual model describing exceedingly slow mineral deposition but it has not been independently confirmed. Assumption 2 -Several equations describing "equilibrium" fractionation in the system calcite-water are available (Craig, 1965;O'Neil et al., 1969;Friedman and O'Neil, 1977;Kim and O'Neil, 1997;Coplen, 2007;Tremaine et al., 2011). The choice of equation affects the calculated temperatures; for example, the application of Coplen (2007) yields temperatures 8-10 • C warmer than those of O'Neil et al. (1969). Assumption 3 -Between ca. 1.0 and 2.7 Ma the isotopic composition of meteoric precipitation in the southern Great Basin was significantly enriched (by 2.5-3.0 ‰) in 18 O relative to modern day values (Winograd et al., 1985;Dublyansky et al., 2011). Meteoric water in the unsaturated zone of Yucca Mountain was similarly enriched in 18 O around 11 Ma (Feng et al., 1999, Sect. 4.1). Indications that δ 18 O of mineral-forming water could have varied in space and/or time are also found in the Yucca Mountain data. Pairing the fluid inclusion temperatures with δ 18 O of host calcite, Wilson et al. (2003, Sect. 7.1) calculated δ 18 O of the parent water to range between −15 and −5 ‰. With similar calculations Whelan et al. (2008, Sect. 5.2) obtained δ 18 O ranging −13 to −7 ‰. Finally, Dublyansky and Spötl (2009) reported that the oxygen in fluid inclusion waters in the Yucca Mountain calcite is significantly enriched (2 to 8 ‰) in 18 O compared to the "normal" meteoric water compositions, apparently due to the waterrock exchange. Because the assumption of constant isotopic properties in waters at Yucca Mountain is not supported by factual data, the merit of using a fixed δ 18 O value for paleotemperature reconstructions for the calcite precipitates there appears to be questionable.
The effect of any or all of the assumptions discussed above being incorrect will be a discrepancy between the calculated and the true temperatures. To check whether this is the case for the Yucca Mountain samples, "isotopic" temperatures calculated using the approach of Whelan et al. (2002) were compared with fluid inclusion temperatures measured from the same calcite (data from Wilson et al., 2003;Whelan et al., 2008). The differences between the measured and the calculated temperatures ranged from −13 • C to +32 • C. In summary, the temporal trends observed in δ 18 O values of calcite do reflect the overall decrease of temperature with time. The large uncertainties associated with the method of calculation, however, render the δ 18 O temperatures unsuitable for quantitative benchmarking purposes.
Best-fit curves -The blue curve in Fig. 2 is the third-order polynomial fit to the data points obtained by isotopic calculations; the other two curves are calculated using the δ 18 O values that are greater and smaller by 2 ‰. The curves, thus, inherit most of the uncertainties of isotopic calculation discussed in previous section. In addition, the curves were constructed in a methodologically erroneous way: the polynomial fit calculation included several data points for which only minimum ages were known (arrows in Fig. 2). The bestfit curves thus have substantial uncertainty for temperatures below 40 • C and become unreliable for temperatures exceeding 40 • C (ages greater than 6 Ma).

Improvement of the benchmark
The currently available data (from Wilson et al., 2003, Fig. 8;Whelan et al., 2008, Table 4) along with associated uncertainties are summarized in Fig. 3. In some cases, inclusions from the same samples and of the same age show fluid temperatures that differ by 15-30 • C (within the dashed oval in Fig. 3). The MICH model describes conductive heating and subsequent cooling of a large rock mass. Abrupt temperature changes are not possible in this model. Because of this conceptual constraint, the observed difference of the temperatures can only be explained by the secondary or pseudosecondary character of the lower-temperature inclusions (i.e., these inclusions were trapped after the formation of the given part of the mineral). Since the time of entrapment for these inclusions is not known, they cannot be used to constrain the thermal history. The endpoint of thermal history (zero time) is constrained by the temperatures expected in the unsaturated zone of Yucca Mountain during different Pleistocene climate states.
The two benchmarks (Figs. 2 and 3) are mutually consistent for times older than ca. 7 Ma. Between 7 and 2 Ma, the benchmark of Whelan and co-authors produces faster cooling (the same temperatures are attained up to 1.5 myr earlier).

Comparison of modeling results with the benchmark
Two versions of the computational modeling of the MICH conceptual model have been reported Whelan et al., 2008). In Fig. 4 the results of this modeling are compared with the improved benchmark. Both versions were run using the HEAT/HEAT3D code (K. Wohletz, © 1998Wohletz, © -2001, The Regents of the University of California). Marshall and Whelan (2001) modeled a 30 km wide, 7 km thick magma body, emplaced at 2.5 km below surface. A partial evaluation of this modeling exercise can be found in BSC (2004a). The conclusion of this study was that modeling failed to match the empirical temperature-time data, which is also obvious from Fig. 4a. The latest results of thermal simulations were reported in Whelan et al. (2008, Sect. 6.3, Fig. 8). The updated model retains most attributes of the earlier model setup, including the geometry and depth of emplacement of the magma chamber. Modifications included (a) hydrothermal convection in the rocks adjacent to the magma chamber, (b) a 100 m thicker unsaturated zone following deposition of the younger Timber Mountain Group tuffs, and (c) subsequent, constant-rate thinning of the overburden due to erosion. It is apparent from the comparison of Fig. 4a and b that the new version of the model brings the temperature-time curves closer to the benchmark, which makes them significantly different from those obtained by Dublyansky and Polyansky (2007, Fig. 10), for a similar model setup and configuration. The authors used the code HYDROTHERM (Hayba and Ingebritsen, 1994), which is better suited for modeling convective two-phase flow through a porous medium than the HEAT3D code used by Whelan et al. (2008).
The technical evaluation documented in Appendix A indicates that the improvements in modeling outcomes were produced by (1) the use of features of the code, which do not render quantitatively realistic solutions (essentially, the code was used beyond the limits of its qualified reliability); (2) arbitrary changes in model parameters; (3) geologically/hydrogeologically unrealistic mesh design; and (4) the usage of non-conservative model parameters. Importantly, all these modifications did not result in agreement between the simulation outcomes and the benchmark; the discrepancies range up to 30-40 • C and 3-5 myr (Fig. 4b).

Evaluation of auxiliary sub-models
Several auxiliary sub-models proposed in the literature in order to bring the computational result of the MICH model into agreement with the benchmark (BSC, 2004a;Houseworth and Hardin, 2010;Whelan et al., 2008) are discussed below. The final sub-model discussed in this section, the cooling action of meteoric infiltration, has not been considered in conjunction with the MICH model before. (4) Lateral distance 7 km, depth ca. 250 m, simulation includes effects of water convection, lateral advective flow of water, and flow in fault zones confining the Yucca Mountain tectonic block (Dublyansky and Polyansky, 2007; HYDROTHERM); (5) Lateral distance 8 km, depth decreases from 300 m at 11.5 Ma to 200 m at present, water convection is allowed (Whelan et al., 2008;HEAT3D). The thick grey line, grey field and circles with error bars show the benchmark as updated in this study.

Continued injection of magma, including injections close to Yucca Mountain
Continued injection of magma into the shallow crust in the vicinity of the Timber Mountain volcanic center after 11 Ma, and/or intrusion of magma closer to Yucca Mountain were proposed in BSC (2004a) to resolve the discrepancy between the computational results of Marshall and Whelan (2001) and the benchmark data. However, geological or geophysical information supporting this hypothesis is lacking (NRC, 2005, p. 16). Before it can be accepted, the purported magma body or bodies would need to be identified, their locations, sizes and time of intrusion would have to be determined, and the thermal effect would have to be numerically modeled. Results compatible with the benchmark so far have been obtained only in simulations where the reference point was located immediately above the margin of the pluton; this means that a hypothetical magma body would have to extend southward as far as Yucca Mountain (the ESF tunnel complex).

Lateral subsurface flow
The lateral subsurface flow of thermal water from the Timber Mountain caldera toward Yucca Mountain was proposed in BSC (2004a) to explain the elevated fluid inclusion temperatures. The unsaturated zone at Yucca Mountain is thought to have formed shortly after ca. 11.6 Ma, and to have persisted since that time (DOE, 2001(DOE, , 2008. The vitric tuffs of the unsaturated zone do not exhibit devitrification or pervasive alteration, as would be expected if heated waters had moved through them en masse for any extended period of time. This means that the "lateral subsurface flow" of thermal waters could only have occurred within the saturated zone; the heating of the unsaturated zone rocks above would still have had to be by conductance. The lateral outflow of thermal water from Timber Mountain is known to have occurred during 10.5 to 11.0 Ma (cf. Sect. 2.2, Fig. 1). Near Yucca Mountain it affected only the deep-seated rocks, below ca. 1000 m. Mineralogical and isotopic evidence (Bish and Aronson, 1993, Sect. "Paleogeothermal conditions"; Feng et al., 1999, Sect. 4.1) as well as thermal modeling (Dublyansky and Polyansky, 2007, Sect. 4.3.3) indicate that this hydrothermal system did not cause heating of the unsaturated zone at the ESF level that is commensurate with the fluid inclusion temperatures.

The presence of additional overburden
The presence of additional overburden, later removed by erosion, would increase the depth and, respectively, the temperatures at the reference points in early stages of the process (BSC, 2004a;Houseworth and Hardin, 2010). It is thought that no more than ca. 100 m of the overburden could have been removed from Yucca Mountain over the last 10 myr (YMP, 1993;DOE, 1998). A 100 m overburden was implicitly included in the simulations of Marshall and Whelan (2001), which calculated the temperatures at a depth of 250 m (BSC, 2004a), whereas the ESF minerals with the highest homogenization temperatures (> 70 • C) were collected from depths of 30 to 80 m. It was then explicitly modeled by Whelan et al. (2008). Sensitivity simulations of Dublyansky and Polyansky (2007) have shown that in order to bring the modeling results into approximate agreement with the empirical data, up to 1000 m of overburden would have to be added. Recently, Houseworth and Hardin (2010) pointed out that the thermal conductivity of the overburden tuffs could have been smaller than that used in sensitivity simulations of Dublyansky and Polyansky (2007) . New simulations, employing this lower thermal conductivity showed that the thickness of overburden still would have to be at least 700 m. The hypothetical addition of such deep overburden would significantly increase the long-term erosion rate; because this is a regulatory concern, acceptance of the "additional overburden" conjecture would require a revision of the Yucca Mountain safety case. Whelan et al. (2008) suggested that vapor-phase convection cells could have developed within the volcanic tuffs above the water table when the water had a high temperature, and speculated that circulation of this vapor heated the unsaturated zone rocks. M. A. Walvoord (US Geological Survey, written communication, 2003) noted that "preliminary simulations of vapor-phase flux in the unsaturated zone above a water table at near-boiling temperatures indicate that vaporphase convection cells could develop within the TSw". Additionally, Whelan et al. (2008Whelan et al. ( , p. 1071) noted that "preliminary simulations indicate that heat transport by vapor-phase convection may account for the high temperatures near the TSw-PTn contact". Because neither the technical details nor the numerical results of these preliminary simulations are provided in the paper, the claim that the proposed mechanism may indeed account for the high temperatures in the unsaturated zone of Yucca Mountain cannot be verified. Conceptually, the proposed mechanism faces at least two serious objections.

Heating by vapor-phase convection
Objection 1 -The mechanism postulates near-boiling temperatures at the water table. There is no geological, geochemical, or mineralogical evidence suggesting that the water table under the repository block has ever been at such high temperatures. During the last large-scale hydrothermal event in the area, the Timber Mountain caldera episode (10.5-11.0 Ma) when the convective outflow occurred under Yucca Mountain, the temperatures at the water table were substantially lower than 100 • C (Bish and Aronson, 1993) (cf. Fig. 1). To make the proposed mechanism acceptable, it must first be reconciled with the hydrogeological history of the area, and independent evidence for near-boiling temperatures at the water table must be presented.
Objection 2 -The elevated temperatures at Yucca Mountain have been determined from fluid inclusions in secondary calcite. This means that the elevated temperatures and the conditions suitable for deposition of calcite must have been present simultaneously. The model describes a convective rise of the hot/warm vapor phase (a mixture of underground air and water vapor) from a hot water table upward into the fractured rock. Because the overlying rock becomes cooler with decreasing distance to the surface, condensation of the water vapor is expected on the fracture walls. The resulting condensate would dissolve the CO 2 present in the subterranean atmosphere and thus become acidic and aggressive with respect to calcite. This mechanism, known as condensation-corrosion, produces sizable chambers in some hypogene carbonate caves (Ford and Williams, 2007, Sect. 7.11;Dublyansky, 2012). According to the general conceptual model of Whelan et al. (2004Whelan et al. ( , 2008, calcite in the unsaturated zone of Yucca Mountain was deposited from small amounts of water infiltrating from the soil zone. Due to the nonlinear relation between the calcium equilibrium concentration and carbon dioxide pressure, the mixing of slightly supersaturated infiltration water with aggressive condensate would necessarily produce water, which is aggressive with respect to calcite (mixing corrosion; e.g., Ford and Williams, 2007, Sect. 3.7.6). Deposition of calcite does not seem to be possible within the proposed vapor-phase convection mechanism.

Elevated heat flows related to extensional tectonics
Houseworth and Hardin (2010) postulated that the thermal history at Yucca Mountain has been affected not only by nearby magmatic activity, but also by the regional heat flux caused by tectonic extension. They contended that the heat flow could have been as large as 300 mW m −2 at 13 Ma and gradually decreased to less than 100 mW m −2 at present, and opined that these high heat flows explain the elevated temperatures in the unsaturated zone at Yucca Mountain. This contention is based on the relationship between the heat flow and the extension rates in the basin and range (Lachenbruch and Sass, 1978, Figs. 9-14) and the history of extension at Yucca Mountain over the last 16 myr (Snow and Wernicke, 2000, Fig. 12). The postulated, extremely high Late Miocene heat flows in the area, however, are in conflict with the sitespecific geological data.
For example, conodonts with a color alteration index (CAI) of 3 were reported from Late Silurian dolostone at a depth of ca. 1800 m in borehole UE25p#1, located ca. 2 km to the east of Yucca Mountain (Carr et al., 1986, App. III). This low CAI indicates that the highest temperature the rock was exposed to since the Late Silurian was 140-180 • C. Assuming that these temperatures were caused by the elevated Late Miocene heat flows, as postulated by Houseworth and Hardin (2010) -although these temperatures could just as well be related to much older, post-Silurian burial of the rock -the corresponding (maximum) heat flows for the area are 110-150 mW m −2 . The paleotemperature gradients that existed during the large-scale hydrothermal event at 10.5-11.0 Ma were reconstructed by Bish and Aronson (1993, Fig. 6) for three Yucca Mountain boreholes from illitesmectite mineralogy and fluid inclusion data (cf. Fig. 1). For the upper ca. 1000 m of the rock mass, the paleogradients are in the 30-36 • C km −1 range, corresponding to heat flows of ca. 50-65 mW m −2 . The hypothesis of Houseworth and Hardin (2010) is therefore not supported by the site-specific data.

Cooling action of meteoric infiltration
Infiltration (and, in the fractured near-surface zone, air movement) leads to an additional cooling of the rocks in the unsaturated zone. The intensity of cooling depends on the infiltration rates; the link between the two parameters is so strong that the temperature measured in boreholes was used to assess the infiltration fluxes at Yucca Mountain (Bodvarsson et al., 2003). The cooling may be slight when the climate is semi-arid and infiltration rates are low, as is the case at Yucca Mountain today (Sass et al., 1988). It should have been substantially greater in the past, when climate in the area was cooler and wetter (Hay et al., 1986;Sharpe, 2007). The cooling action of infiltration can create a rain curtain state of near-isothermal conditions from the surface to depths of 800-1000 m. Such an effect has been observed at one modern geothermal system in Oregon (Swanberg and Combs, 1986). A rain curtain occurred in the Yucca Mountain unsaturated zone during the Timber Mountain caldera hydrothermal episode 10.5-11.0 Ma (Bish and Aronson, 1993, Sect. "Paleohydrologic conditions"). The inclusion of the effect of cooling by infiltration would lower the temperatures computed for the MICH model in the depth interval of interest (100-300 m). In fact, the temperatures in this zone could have been entirely controlled by this effect, as illustrated in Fig. 1. 6 Additional approaches to validation of the MICH model

The analog-system observations
Observations from natural analogs are an integral part of the model development process in programs related to geological disposal of nuclear waste (DOE, 2004). The term "natural analog" refers to natural systems in which processes similar to those expected to occur in a nuclear waste repository are thought to have occurred. Analog investigations may determine the conditions under which the processes occur and the effects of the processes, as well as the magnitude and duration of the processes (Simmons, 2003). Natural analog studies represent a somewhat weaker form of validation: while they may confirm that an hypothesized process is generally possible and has occurred elsewhere, this does not necessarily imply that the process has occurred at the site which is being evaluated. The Long Valley caldera in California has been proposed as a useful analog for examining the processes postulated in the MICH model (BSC, 2004a). The parameters of the parent magma chamber at Long Valley appear to be similar to those of Timber Mountain (an elliptical form measuring 20 km × 30 km and 2 km thick, with a depth of emplacement of 5-7 km (Bailey et al., 1976); the estimated magma temperature is ca. 800 • C (Hildreth and Spera, 1974). Silicic volcanism at the Long Valley igneous system began ca. 2 Ma and continued intermittently until the time of the caldera collapse at 0.7 Ma (Bailey et al., 1976). Shortly after the collapse a resurgent dome formed, and silicic and intermediate volcanic material has been discharging around it virtually up to the present time. Heat flow was measured in 11 shallow boreholes (150-300 m) in the vicinity of the caldera by Lachenbruch et al. (1976), who concluded that although the magma chamber has likely been replenished during its 2 myr long eruptive history, no conspicuous indications of thermal anomalies related to this chamber could be detected outside the caldera margins. This conclusion was upheld in BSC (2004a): ". . . a thermal anomaly outside of the Long Valley caldera would not be detectable 700 000 years after the last major phase of magmatic activity" (p. H-21). Observations at this natural analog system, therefore, do not support the validity of the MICH model, which describes a strong thermal anomaly at a distance of several kilometers from caldera margin existing for several million years.

Observations on spatial structure of the thermal field
In previous sections we evaluated the temperature-time relationships, while the spatial information was reduced to a single parameter, distance from a hypothetical magma body. Another potentially relevant feature is the spatial structure of the thermal field. The field expected in the MICH model is characterized by a strong lateral gradient with the highest temperatures immediately adjoining the magma body and decreasing rapidly with increasing distance away from its margins. A critical validation test would consist of a comparison of the actual spatial distributions of paleotemperatures around the proposed heat source, the Timber Mountain magma body, and those expected from the model. The temperature data obtained from the ESF tunnel cover an area of a few square kilometers which, in the evaluations of the MICH model given above, has been considered as a single point for the sake of simplicity. The distribution of maximum fluid inclusion temperatures within this area exhibits a pronounced east-west gradient, orthogonal to the gradient expected within the MICH model (Fig. 5). This distribution requires that the source of heat be located to the east of the repository block.
Fluid inclusion homogenization temperatures from calcite from borehole USW G-2 located ca. 3 km to the north of the ESF complex (i.e., closer to the presumed heat source, where higher temperatures would be expected) are reported in USGS (2002) and Whelan et al. (2008, Wilson et al., 2003;Whelan et al., 2003Whelan et al., , 2008. have been farther to the east than previously thought". Whelan et al. (2008) did not discuss the significance of these data for the MICH model. The idea of "relocating" the magmatic heat source to the east of Yucca Mountain is in conflict with the site-specific geological evidence because no significant buried magma bodies are known in this direction. The low conodont color alteration index of the Late Silurian rocks penetrated by the UE25 p#1 borehole ca. 2 km to the east of Yucca Mountain (Carr et al., 1986) speaks strongly against such a possibility. The spatial structure of the paleotemperature field, as recorded by fluid inclusions, effectively falsifies the MICH model: the temperatures do not increase toward the presumed heat source, Timber Mountain, as they should. (For discussion of alternative interpretations put forth to explain the paleothermal gradients in the unsaturated zone of Yucca Mountain see Appendix B).

Discussion and conclusions
The development and validation of scientific models associated with geological disposal of nuclear waste in the United States are regulated by the US Department of Energy's Quality Assurance Requirements and Description (DOE, 2004) and by a set of administrative procedures (e.g., AP-3.10Q -Analyses and Models, AP-3.11Q -Technical Reports, and AP-SI.1Q -Software Management). The analysis presented in this paper demonstrates that none of the three approaches mandated by regulations to validate the model, computational modeling, analog-system observations, and observations on the structure of the paleotemperature distribution support the MICH model. The model, therefore, cannot be considered validated. Nevertheless, the model has become firmly established in the scientific record. This was aided by the (erroneous, as argued in this paper) notion that the conceptual model of the conductive heating at Yucca Mountain has been validated by means of computational modeling. Statements to this effect have appeared in publications by the authors of the model Whelan, 2000, p. A-259, 2001, p. A-375;Whelan et al., 2002, Sect. 5.3;2003, Sect. IV;2004, Sect. 3.2, 2008Marshall et al., 2005, Sect. 2;Paces and Whelan, 2012, p. 245; see also Appendix C). Other researchers relied on the MICH model to interpret their data (Wilson et al., 2003, Sect. 7.6;Bryan et al., 2009, Sect. 3). The model was incorporated, implicitly or explicitly, into the DOE program documents, such as the Yucca Mountain Science and Engineering Report (DOE, 2001, p. 4-402) and the Yucca Mountain Site Description (BSC, 2004b, p. 7-81). It served as one of the key arguments for the exclusion of the FEP (feature, event, process) hydrothermal activity from consideration in the Yucca Mountain Total System Performance Assessment (BSC, 2004c;SNL, 2008;cf. Dublyansky, 2007).
Fluid inclusions in secondary minerals in the unsaturatedzone of Yucca Mountain provide direct evidence that thermal waters accessed this zone in the past. The flow of thermal water is a process that could compromise the performance and the safety of the nuclear waste disposal facility. The MICH model is the only conceptual model considered in the DOE Yucca Mountain technical documentation to explain the presence of thermal waters in the unsaturated zone of Yucca Mountain. The effort to validate the model documented in this paper was unsuccessful. As a result, the DOE's Yucca Mountain safety case does not have plausible explanation for this natural phenomenon, potentially detrimental to the safety of nuclear waste disposal. It is important, therefore, that the reliance on this model is discontinued, and the appropriateness of decisions based on this model is re-evaluated.

A1 Introduction
The initial version of the MICH model (designated here MICH_00, Table A1; Marshall and Whelan, 2000) was not described in sufficient detail. A subsequent version of the model (MICH_01), for which numeric simulations were performed, included conduction of heat from shallow crustal pluton and water convection, which was limited to the area located above the magma chamber. This version of the model was evaluated by means of numeric modeling (Marshal and BSC, 2004;Dobson, 2003), whose results did not match the benchmark (a composite temperature-time curve, based on empirical data -fluid inclusion temperatures and radiometric age dates for secondary minerals at Yucca Mountain). More sophisticated simulations, which included effects of advective and convective water movement in the model domain (Dublyansky and Polyansky, 2007), were equally unsuccessful (Fig. 4a).
The latest version of the model (MICH_02; Whelan et al., 2008) purportedly included heat transfer by water convection in rocks surrounding the magma chamber. In contrast to previous version, simulations produced "cooling histories" plotting relatively close to (although still not matching) the benchmark (Fig. 4b).
In this Appendix I evaluate the technical aspects of modeling performed by Whelan and co-authors on MICH_02 in order to assess whether or not the results can be accepted as viable. The pertinent questions are: (1) Can the numeric modeling results be reproduced? (2) What changes in model setup and parameters and the usage of the code produced such a dramatic change in the modeling results? and (3) Are the implemented changes reasonable/acceptable? This evaluation was initially hampered by lack of information. The original publication of Whelan et al. (2008) does not contain sufficient technical data about model and modeling process. The initial request for information made to the authors was declined. Because of that, the evaluation was not included in the first submission of this paper. The necessary materials (Supplement 1) were eventually obtained with the help of Elsevier Publishers and the Office of Science Quality and Integrity of the US Geological Survey.
For simulations of MICH_02, Whelan et al. (2008) used the code HEAT3D, ver. 4.10.0517. The replication and sensitivity simulations discussed in this Appendix were performed using the more recent ver. 4.11.0533 of the code.

A2 Reproducibility simulations
By using the mesh design, model parameters, and scenarios of magma emplacement and permeability changes of Whelan et al. (2008), provided in the Supplement 1, the "cooling curves" were simulated for the reference point, located 8 km from the boundary of the model pluton, at a depth of 300 m. The obtained cooling curve had a stable shape, which was reproduced in multiple realizations, run on different computers, with the mesh created anew before each simulation. The curve obtained has an overall shape similar to that presented in Whelan et al. (2008) but with different parameters (Fig. A1). The modeled temperature was up to 18 • C lower (e.g., between 11 and 10 Ma) and the rock cooled down to the same temperature (e.g., 45-50 • C) up to 1.9 million years earlier.
The major difference pertains to the height of the "plateaus" on the descending limb of the cooling curve, which were lower in the replication simulations. The heights of these plateaus (e.g., 78, 50, and 30 • C) are controlled by the porosities assigned to the rocks (the higher the porosity, the higher the temperature). The abrupt drops following the plateaus (e.g., at 10 and 8 Ma) reflect decreases in porosity prescribed in the model scenario (cf. Sect. A4.1). I was unable to reproduce the plateaus of Whelan et al. (2008) by increasing the porosity up to f = 0.4, the maximum value allowed in HEAT3D. The replication modeling was thus only partly successful. The reasons for the discrepancy between the cooling curves obtained by Whelan et al. (2008) and in this study remain unclear.
In subsequent sections I will evaluate what controls the shape and the parameters of the cooling curves shown in Fig. A1.

A3 Analysis of mesh design
The mesh used by Whelan et al. (2008) is shown in Fig. A2. The model domain is capped by a thin (0.5 km) impermeable layer representing the unsaturated zone (UZ) rocks. The heat transfer in this layer is by conduction only. Water convection is allowed in the underlying mesh domain designated "volcanic convecting". Parameters of convection in this domain are controlled by model scenario: the porosities decrease stepwise, at specified times, from f = 0.1 to 0.05 to 0.02 to 0. Heat transfer in all other rocks is by conduction. Magma bodies are emplaced twice at a depth (to the roof) of 2.5 km. The reference point is located at x = 7 km, corresponding to lateral distance from the pluton of 8 km.
One puzzling feature is the presence of an impermeable block ("volcanic non-convecting") in the upper-left part of the mesh. This "hydraulic barrier" does not seem to be connected to the "real-world" saturated-zone hydrogeology of Yucca Mountain. At least no such rocks appear in any of the site-scale hydrogeological models of the saturated zone of Yucca Mountain (e.g., Zyvoloski et al., 2003;Belcher et al., 2012). The effect of this impermeable block on the simulation outcome will be evaluated in Sect. A4.2.  (2000) MICH_01 A 50 km × 30 km mesh with 0.5 km × 0.5 km grid elements. Layer-cake stratigraphic model (from the surface down: volcanic, carbonate, mafic igneous, and metamorphic rocks). 30 km diameter disc-shaped pluton (volume 5000 km 3 ; T = 900 • C) emplaced at depths of 5 or 2.5 km. Hydrothermal convection directly above the magma chamber is allowed. The cooling histories (T vs. time) were simulated for reference points located 0 and 4 km away from edge of the pluton; depth 250 m.  In physical terms, the model shown in Fig. A2 is similar to a pressure cooker. A body of water-saturated permeable rock is confined by impermeable "walls" and a "lid", and is heated from below by a "heater" -geothermal gradient. Two additional heaters (plutons) are introduced into the "vessel" during the simulation. Geometrical setting of the model is ideal for the development of the Rayleigh-Bénard convection, which would operate as long as the heat supply continues. From the standpoint of hydrogeology, however, the model is not realistic, as it does not consider the advective recharge and discharge. Both processes would be expected to speed up the cooling of the system because the recharge introduces cool water and the discharge removes heated water from the system.

A4.1 The effect of water convection
When water convection is enabled in simulations, the temperature at the reference point quickly increases and then remains constant (creating a "plateau" on the simulated temperature curve) as long as the specified porosity value remains unchanged. Sensitivity simulations show that relatively high temperatures (> 60 • C) can be produced at the reference point by a normal geothermal gradient of 30 • C km −1 , if sufficiently high values of porosity are assumed (Fig. A3).
The question remains: are the intensities of convective heat transfer and the resulting high temperatures at the reference point realistic? To answer this question, it is important to understand the limitations of the code. HEAT3D provides very simplistic convection solutions, which have not been quantitatively validated. The code does not model the Table A2. Comparison of parameters used in thermal simulations of the MICH model by Marshal and Whelan (2001) and Whelan et al. (2008) with realistic parameters for the Timber Mountain volcanic center.
Initial temperature of magma 900-1000 • C 700-750 • C (upper part) to 900-950 • C (lower part) (Mills et al., 1997;Bindeman and Valley, 2003) Residence time of magma Permanent emplacement(no evacuation) Between 10 000 and 100 000 years for early plutons, followed by eruptive evacuation (Bindeman and Valley, 2003). Permanent for the post-11.45 million-year-old resurgent dome. * Structural boundaries of calderas provide the best estimate of the shapes of the underlying magma chambers (e.g., Smith and Shaw, 1978). Whelan et al. (2008, Sect. 6.3) speculated that the pluton under the Timber Mountain caldera could have had a subsurface footprint that extended beyond the caldera margins (i.e., closer distance to reference point). There is no site-specific information supporting this conjecture. Rayleigh-Bénard convection; instead it approximates Darcy conditions of porosity-dependent flow in permeable rock that carries heat down the thermal gradient. Due to these constraints of the code, the movement of fluids (i.e., hydrogeological situation), and the associated heat transfer, were not (and could not have been) realistically reproduced in the simulations of Whelan et al. (2008).

A4.2 The effect of the "intrusion advection" option
In the next series of the sensitivity simulations, the emplacement of plutons was modeled with water convection disabled. The results of the conduction-only simulations strongly depended on whether or not the option intrusion advection was enabled (Fig. A4). With this option turned off, two consecutive emplacements of magma bodies produced very minor heating at the reference point (curve 2 in Fig. A4). In contrast, when "intrusion advection" was enabled, the emplacements produced two abrupt, peak-shaped temperature increases at the reference point, with the second peak being almost twice as high as the first one (curve 1 in Fig. A4). Considering that in this simulation the heat transfer is supposed to occur only by conductance, the sharpness of the peaks is puzzling. Meanwhile, from our reproducibility simulation it is clear that the high peaks distinguishing the new cooling curve of Whelan et al. (2008) are the result of the use of this option.
To assess the plausibility of the results, it is instructive to examine the function of the intrusion advection option in HEAT3D. This option models the effective heat transport caused by displacement of surrounding rock by the intrusion of magma. Adding a new magma body instantaneously (i.e., during addition of new magma into the mesh) transports the preexisting heat (temperature) outward and upward in the mesh. The location to which this preexisting rock is moved depends upon the shape of the new intrusion; a quasi equidimensional intrusion, such as the one used in this model, pushes rock both upward and outward. The code redistributes the pre-existing heat into a volume of rock around the intrusion, assuming some mixing with cooler ambient rocks.
Because the HEAT3D documentation does not explain how exactly the mixing is performed, this was evaluated through simulations. Modification of the thermal field at the moment of pluton emplacement is shown in Fig. A5a. If the intrusion advection is disabled, the temperature in the rock surrounding the pluton is controlled by the geothermal gradient. In contrast, when the intrusion advection is enabled, x (km) Figure A2. Mesh used by Whelan et al. (2008). The temperatures were modeled for the reference points, located at z = 200-300 m from surface (lateral position is marked by a flag; not to scale). Water convection was allowed in the "volcanic convecting" domain.  Figure A3. Sensitivity simulation illustrating the role of water convection. The graph shows the temperature at reference point (z = 300 m). Porosity changes with time as in simulation of Whelan et al. (2008). There was no magma emplacement. the temperature at the pluton boundary and above is increased, by ca. 100 • C. In a horizontal direction, the temperature decreases linearly from the pluton boundary, reaching the "geothermal" values at the boundary of the mesh. At z = 1 km under the reference point (x = 7 km) the instantaneous increase of the temperature is ca. 50 • C, resulting in a geothermal gradient of ca. 80 • C km −1 . If water convection is enabled, this instantaneous effect of the intrusion advection is further augmented (Fig. A5b).
The question of to what extent the intrusion advection option, as implemented in HEAT3D, reflects the real-world thermal and rheological effects of pluton emplacement is non-trivial. In modeling heat flow around large magma bodies, a fundamental consideration must be made. Can a model where a large magma body appears instantaneously within host rocks be realistic? What happens to the rock that the magma displaces? For realistic calculations, one must somehow model the fact that the magma body takes a long time to grow, a time during which heat is supplied to the host rock and host rock is displaced.
Because most of the improvement of the match between the new cooling curve of Whelan et al. (2008) and the benchmark is due to the intrusion advection option, it is critically important that the potential uncertainties associated with this decision are discussed. This has not been done in the original publication by Whelan et al. (2008). Actually, the very fact that this option was enabled in thermal calculations was neither mentioned in the paper nor documented in the materials provided for this evaluation by the authors of the model (see Supplement 1). (2) the intrusion advection is disabled. There was no water convection.

A4.3 The effect of the impermeable block
As was noted above, the geological/hydrogeological nature of the block of impermeable rock, confining the convective domain at the upper-left margin of the mesh (Fig. A2) is not clear. In order to assess its role, sensitivity simulations were run with and without this "hydrological barrier". Simulations show that the inclusion of the barrier increases the temperature reached at the reference point to between 10 • C (at z = 100 m) and 32 • C (at z = 300 m). The effect of the impermeable rock on augmenting the cumulative effect of the intrusion advection and water convection is also apparent in Fig. A5b.

A4.4 The thickness of the unsaturated zone
The model setup of Whelan et al. (2008) includes the UZ layer, which is only 500 m thick. The modern day thickness of the UZ in the Yucca Mountain area is between 500 and 700 m (Montazer and Wilson, 1984;Wu et al., 2002). In addition, some 100 m of the rock could have been eroded away during the 12 million-year-long history of the mountain (this was also accounted for by Whelan et al., 2008). Sensitivity simulations show that with increasing thickness of the UZ, the temperatures at the reference point become lower and less responsive to the convective heat transfer in the underlying model domain (Fig. A6). Horizontal temperature profiles at z = 3000 m, capturing the consecutive events in a simulation following the scenario from Supplement 1: (5) conductive heat transfer; (6) water convection allowed; (7) emplacement of the pluton #1; (8) effect of water convection 50 ka after pluton #1 emplacement; (9) emplacement of the pluton #2; and (10) effect of water convection 50 ka after pluton #2 emplacement. The vertical dashed line corresponds to the x coordinate of the reference point.

A4.5 Conservativeness of modeling
The question remains as to whether the mismatch between the computational outcomes of Whelan et al. (2008) and the benchmark can be reduced by using less conservative model parameters. Analysis shows that this is very unlikely, primarily because the parameters used for computing the thermal histories shown in Fig. 4b were non-conservative already. One example of non-conservativeness is given in Sect. A4.4. The major source of uncertainty in the MICH model, however, pertains to poorly known parameters of the plutons (dimensions, locations, depths of emplacement, residence time, etc.). Estimates of these parameters for the Timber Mountain caldera system depend on the petrogenetic model selected, which varies from that of a single magma chamber (Lipman et al., 1966), to successive emplacement of magma chambers , to multiple magma batches (Schuraytz et al., 1989;Huysken et al., 1994;Cambray et al., 1995;Mills et al., 1997), to vertically separated magma bodies and re-intrusion (Bindeman and Valley, 2003). Temperatures, emplacement dynamics, depth, and residence times of magma will all vary dramatically between these models. For example, a very shallow depth of emplacement of 2.5 km was used in the simulations shown in Fig. 4; greater emplacement depths decrease the simulated temperatures at the reference points. However, although the depth of 2-3 km was proposed for the Timber Mountain plutons on the basis of geometrical consideration of the calderas (Byers et al., 1976), these estimates were subsequently revised by Warren et al. (1989) and Mills et al. (1997) who, using mineral geobarometers, estimated the depth of emplacement as 10 km to 12.5 km. Comparison of selected parameters in Table A2 shows that the modeling of Marshal and  and Whelan et al. (2008) was non-conservative and likely overestimated both the temperatures and the cooling times.

A5 Conclusions
This evaluation demonstrates that the improved match between the cooling curves of Whelan et al. (2008) and the benchmark was produced by (1) the use of the computer code beyond the limits of its qualified reliability (convective heat transfer, Sect. A4.1); (2) the use the code features which do not provide quantitatively qualified solutions (intrusion advection, Sect. A4.2); (3) setting model parameters arbitrarily (schedule of porosity changes specified in model scenario, Sect. A4.1); (4) mesh design which is hydrogeologically unrealistic (the presence of the impermeable barrier, Sect. A4.3); and (5) the non-conservative model parameters (Sects. A4.4 and A4.5; see also Dublyansky and Polyansky, 2007). Each modification of the model implemented by Whelan et al. (2008) shifted the computed results either toward higher temperatures or toward longer cooling times, i.e., toward the benchmark.
HEAT3D is a simple and robust research tool, useful for the first-order modeling of relatively simple geological systems. For more complex situations, particularly those involving hydrogeological considerations the code either provides approximate solutions, which are not quantitatively validated, or do not provide solutions at all. For such complex systems, HEAT3D should not be used. The code has not been qualified for use in any DOE's Yucca Mountain Project research, subject to Quality Assurance Requirements and Description (DOE, 2004).
Geological, hydrogeological, and geometric settings of the MICH conceptual model are complex. To be considered realistic, the thermal simulations, among other things, must consider: (1) the conducive heat transfer in the unsaturated zone; (2) conductive, convective and advective (lateral) heat transfer in the underlying saturated zone; (3) impact of the emplaced magma chambers on the thermohydrology in the surrounding rocks; (4) topography of the Earth's surface; etc. While some of these processes can be modeled by HEAT3D with a reasonably high degree of fidelity (e.g., 1), others can only be captured through more or less crude approximations (2 and 3), and some cannot be modeled at all (4). Simulation results reported by Whelan et al. (2008), therefore, cannot be viewed as "realistic" in a quantitative sense; rather, they document the outcome of a computational experiment which has few ties to observed reality.

Appendix B: Evaluation of the paleothermal field at Yucca Mountain
When evaluating information regarding the paleothermal field in the Yucca Mountain unsaturated zone, one needs to address two questions: (1) are there indications of spatial gradients in the paleothermometric data? And (2) if yes, how may these gradients be interpreted? 1598 Y. V. Dublyansky: Evaluation of the US DOE's conceptual model of hydrothermal activity

B1 Observation of spatial paleotemperature gradients
The non-uniform character of fluid inclusion paleotemperatures in samples from the ESF tunnel was first noted by Dublyansky (1998): In terms of the spatial distribution of measured homogenization temperatures, the following observation may be important. Two samples that yielded temperatures higher than other samples. . . are both from Tiva Canyon tuff. Also, both these samples are from the eastern part of the exploratory block. . . closest to the Paintbrush (∼ 2 km to the east of the repository block) fault zone which might have served as major avenue for upwelling fluids. (Sect. 6.8) The observation was based on the limited number data available at that time and was, therefore, considered tentative: "Although it is premature to make strong conclusions on the basis of only two samples, this hypothesis needs to be addressed in the future, when the spatial structure of the ancient upwelling system is studied." (Dublyansky, 1998, Sect. 6.8) More data became available in the course of the University of Nevada, Las Vegas Yucca Mountain Thermochronology project in 1999-2000. In the final report from this Project, Wilson and co-authors asserted that "temperatures recorded across the Yucca Mountain repository horizon do not exhibit a central hot plume and large lateral thermal gradients that are present in geothermal and epithermal systems . . . The lack of a significant temperature gradient and presence, instead, of relatively uniform temperatures argues against an upwelling hot fluid model". (Wilson et al., 2002, Sect. 6.1).
In the follow-up publication, a similar statement appears (Wilson et al., 2003, Sect. 7.3). These statements purport that the early tentative information on the non-uniform character of the paleothermal field at Yucca Mountain was not confirmed. Further, Wilson et al. (2003) proposed the presence of "spatially localized high-temperature fluids" near the north portal of the ESF; no interpretation was offered as to the nature of these localized fluids. In contrast, gradients in paleothermometric data were reported by Whelan et al. (2003; Sect. III.A): "fluid inclusion T h [temperatures of homogenization] in calcite decrease from east to west along the ESF north ramp in the north bend (from about 90 • C to about 60 • C), from northeast to southwest along the ECRB Cross Drift from about 60 to 50 • C, and from north to south along the ESF main drift from about 60 • C to about 45 • C".
One reconstruction of the paleotemperature field based on fluid inclusion data is shown in Fig. B1. It must be emphasized that homogenization temperatures shown on this graph were obtained from the basal, oldest parts of mineral crusts and, therefore, characterize the early stages of mineral deposition. The temperatures decreased with time and eventually attained the ambient values.
Additionally, an important observation is that the distribution of δ 18 O values in early calcite in the repository block virtually mimics that of the paleotemperatures (Fig. A2). This can be explained by the temperature-dependent isotopic fractionation in the CaCO 3 -H 2 O system: "the fractionation of oxygen isotopes between calcite and water is a function of temperature, with the calcite δ 18 O value increasing as temperature decreases. The trends in the δ 18 O values of earlystage calcite. . . are, therefore, consistent with the temperature trends displayed by the FIA [Fluid Inclusion Assemblage] T h measurements" (Whelan et al., 2003, Sect. III.B) In summary, although one may put forth different interpretations, and attach different significance to the non-uniform thermal field revealed by fluid inclusions and stable isotopes, characterization of the data as exhibiting "lack of a significant temperature gradient and presence, instead, of relatively uniform temperatures" (Wilson et al., 2002, Sect. 6.1) is clearly unwarranted. Dublyansky and Smirnov (2003, Sect. 2.4) and Dublyansky et al. (2005, Sect. 2) considered the temperature gradients apparent in Fig. B1 as a true feature of the paleothermal field, characteristic of the early stages of secondary mineral deposition in the Yucca Mountain unsaturated zone. Alternative interpretations are discussed in subsequent sections.

B2.1 Incomplete sampling
The idea that the gradients are apparent rather than real, and result from the failure to find samples with high-temperature inclusions in the part of ESF where only low homogenization temperatures were measured, was informally put forth in a number of discussions. Samples for the fluid inclusion studies were collected from the 7.8 km long ESF tunnel by three research groups in the course of several sampling campaigns. Most of the ESF sites bearing secondary mineralization have been sampled (cf. Wilson et al., 2003, Fig. 2;Whelan et al., 2008, Fig. 2a). The currently available fluid inclusion database consists of ca. 7000 measurements obtained from more than 130 sampling sites (Whelan et al. 2008, Sect. 5.2). The δ 18 O database exceeds 2000 measurements. Independently obtained fluid inclusion temperatures and stable isotope data are remarkably consistent (Fig. B2). The possibility that some occurrences containing either the high-temperature fluid inclusions or characteristically low δ 18 O values were not sampled in the north-south drift of the ESF (where the paleothermal field shows the lowermost values) appears therefore remote.

B2.2 Temporal rather than spatial gradients
With respect to the non-uniform temperatures of early fluids,  opined that:  Wilson et al. (2003), and Whelan et al. (2008). Red crosses in (c) indicate temperatures older than 8.5 Ma (data from Wilson et al., 2003;Whelan et al., 2008). Each mode shown in (c) was calculated on the basis of tens to hundreds of individual measurements reported in , Wilson et al. (2003) and Whelan et al. (2003). Vertical dashed lines correspond to the bends in the ESF (cf. Fig. 5). This hypothesis is supported only by tenuous argumentation. At Yucca Mountain, secondary minerals are found in two types of openings in the rock: lithophysal cavities and fractures/breccias. Lithophysal cavities formed shortly after the emplacement of host tuffs and were subject to vaporphase alteration; the surfaces of some early fractures are also altered. Fractures which bear no trace of vapor-phase alteration should have formed later than lithophysal cavities. Wilson and Cline (2002, Sect. 7.2) stretch this simple logical premise to argue that "secondary minerals in fractures and breccias began to precipitate later than secondary minerals in lithophysal cavities".
The previous statement is a non sequitur. The vapor-phase alteration ended shortly after the emplacement of the ash material, its compaction and conversion into welded tuffs. The complete cooling of the 350 m thick sequence of Topopah Spring ash-fall tuff from its estimated temperature of emplacement to ambient temperatures took about 7000 years (Buesch and Riehle, 2007). The vapor-phase alteration took only a fraction of this time. A fracture which was not affected by vapor-phase alteration could be just a few thousand years younger than a cavity affected by it; a tiny fraction of their more than 12 myr long history.
The contention of  is also not supported by the data. The U-Pb dating results (Neymark et al., 2002;Wilson et al. , 2003Whelan et al., 2008) indicate that (a) all secondary minerals at Yucca Mountain postdate the vapor-phase alteration by 1 to 1.5 myr and (b) many fractures contain minerals, which are as old as their counterparts from the lithophysal cavities (Fig. B2b). Figure B2b also shows that areas in the ESF characterized by different early paleofluid temperatures were accessed by these fluids at essentially the same time (within the error of the U-Pb dating method).
The paragenetic age of calcite at Yucca Mountain can be assessed on the basis of its stable isotope properties. Carbon is a conservative component, whose isotopic composition reflects that of the source(s) of dissolved bicarbonate. Fractionation of carbon isotopes is affected little by temperature (e.g., the fractionation coefficient for HCO 3(aq) -CaCO 3 changes by 0.4 ‰ between 20 to 90 • C; Deines et al., 1974). In the Yucca Mountain samples δ 13 C decreases from 8 to 10 ‰ in the earliest calcite to −8 to −10 ‰ in the latest one (Wilson and Cline, 2002, Sect. 6.4;Whelan et al., 2002). The earliest calcite with characteristic, strongly positive δ 13 C values is present throughout the repository zone (Fig. B2e), from areas near the north and south portals where the highest paleotemperatures were measured to the N-S drift. This is consistent with the conclusion of Whelan et al. (2002, Sect. 4.4): "The large range of δ 13 C values, as plotted against location in the ESF . . . , shows that the entire paragenetic sequence is present in mineral coatings throughout the ESF". Whelan et al. (2003, Sect. IV) and Whelan et al. (2008, Sect. 6.2) suggested that the two highest temperature samples from the first 500 m of the ESF tunnel, one of calcite and one of fluorite, were formed as a result of the transient fumarolic activity. The authors noted the presence of the thin bleached rims on the fractures hosting the minerals at these two locations and interpret these rims as an indication that the fractures transported hot, fumarolic fluids produced during the initial cooling of the Tiva Canyon tuff. Although the evidence presented in support of the fumarolic origin of these samples is by no means definitive, we note that even if the highest temperature calcite sample is excluded (fluorite thermometric data were not used for paleothermometric reconstructions because of their lower reliability at Yucca Mountain), the east-west gradient remains (Fig. B2c). The possible fumarolic activity is thus not relevant for the gradients of the paleothermal field in the Yucca Mountain unsaturated zone.

B2.4 Link with zeolite alteration
Describing the paleotemperature gradients observed in their fluid inclusion and δ 18 O data Whelan et al. (2008) noted: "similar temperature trends are suggested by the intensity of zeolitization in the upper Calico Hills Formation (Bish et al., 2003, Fig. 8). . . . this temperature trend more likely reflects a lateral, east-west, temperature gradient, perhaps correlative with zeolitization intensity in the Calico Hills Formation . . . " (Sect. 5.2.1).
Genetic implications of the apparent correlation between the temperatures measured in secondary minerals in the ESF and the abundance of zeolites in the layer of vitric tuffs located several hundred meters deeper are not discussed by Whelan and co-authors. Without such a discussion, the above statement represents a causal fallacy (correlation does not imply causation).
It is to be noted that the zeolitization of the Calico Hill formation likely occurred between 12.9 Ma (the age of the tuff; Sawier et al., 1994) and11.3 Ma (Broxton et al., 1987, p. 101). As was stated above, none of the secondary minerals from ESF yielded a U-Pb age exceeding 10 Ma, implying a more than 1 myr long gap between the two processes. Further, it was demonstrated that the latest large-scale hydrothermal system (10.5-11.0 Ma; Bish and Aronson, 1993), which likely produced zeolitization, only affected the deep parts of the Yucca Mountain rock sequence (> 1000 m). During this hydrothermal event, the temperatures in the unsaturated zone remained below those indicated by fluid inclusions (Bish and Aronson, 1993;Feng et al., 1999, Sect. 4.1;Dublyansky and Polyansky, 2007, Sect. 4.3.3).

B3 Conclusions
Paleothermal gradients revealed by fluid inclusion homogenization temperatures and the δ 18 O compositions of early calcite must be considered real unless and until convincing evidence otherwise is presented.
No coherent explanation of the origin of gradients can be found in the literature (cf. "localized high-temperature fluids" of Wilson et al. (2003), "fumarolic activity" of Whelan et al. (2003), and trends "correlative with zeolitization intensity in the Calico Hills Formation" of Whelan et al., 2008).
The observed east-west gradients in the ESF area require the heat source to be located to the east of the ESF. This is in direct conflict with the MICH model, in which the source of heat is located to the north of the ESF, under the Timber Mountain caldera. The observed structure of the paleothermal field is a strong argument against the MICH model.
An argument can be made that the exact configuration of the paleothermal field depicted in Fig. B1 is uncertain because for some areas in the ESF it is constrained by a limited number of measurements; therefore, different configurations of isotherms can be drawn (Wilson and Cline, 2005, Sect. 2.1). This argument is legitimate. The paleothermal field shown in Fig. B1 must be viewed as a first-order feature only. This, however, does not diminish its significance.
Any further development of this line of enquiry would have to include fluid inclusion and stable isotope studies of secondary minerals from numerous boreholes drilled within and outside the ESF footprint. In addition to providing better lateral coverage, such studies would provide sorely needed information on distribution of temperatures with depth.

C1 Introduction and promotion of the MICH model
The meteoric infiltration -conductive heating (MICH) thermal model was introduced by the USGS researchers in 2000 and 2001 at the annual meeting of the Geological Society of America. Two short abstracts Whelan, 2000, 2001) mention thermal modeling and state that the results of simulations are in agreement with the age, and empirical thermometric data obtained from secondary minerals at Yucca Mountain (emphasis is added in all quotations below): This trend indicates a gradual cooling of the rocks over millions of years, in agreement with thermal modeling of magma beneath the 12-Ma Timber Mountain caldera just north of Yucca Mountain. This model predicts that temperatures significantly exceeding current geotherm values occurred prior to 6 Ma." (Marshall and Whelan, 2000, p. A-259) The simulations indicate that modern geothermal gradients were reached at 6 Ma to 3 Ma. These results are in general agreement with paleotemperature data from fluid inclusions and isotopic compositions of secondary calcite at Yucca Mountain." (Marshall and Whelan, 2001, p. A-375) On 9 May 2001, J. Whelan testified at the US Nuclear Waste Technical Review Board meeting in Arlington, VA. In his presentation, he again stated that thermal simulations were successful: "so, to conclude, both fluid inclusions and calcite delta O-18 indicate elevated temperatures during the early and intermediate stages of calcite formation. Those temperatures are consistent with a likely thermal history of the unsaturated zone tuffs as indicated by the age constraint temperature data and by thermal modeling." (NWTRB, 2001, p. 151) The record demonstrates, however, that contrary to these repeated statements, simulations carried out in 2001 failed to reproduce the empirical temperatures and times. The document (BSC, 2004a), in which the USGS modeling effort of 2001 was evaluated states: Marshall and Whelan (2001) concluded that the presence of a long-lived magma chamber at the Timber Mountain volcanic center could account for elevated thermal conditions in the vicinity of the repository up to around 6 Ma.
However, closer evaluation of the model results . . . indicate that the magmatic activity at Timber Mountain as represented by these simulations would only produce minor and relatively shortlived thermal perturbations for the Yucca Mountain area. (p. H-9-H-10).
Despite the failure of the simulations to support their model, between 2001 and 2012, the authors of the MICH model continued to make statements regarding the success of thermal modeling: Warmer depositional temperatures in the past reflect the prolonged thermal input to the UZ from the ongoing regional magmatic activity . . . Yucca Mountain tuffs were erupted between 15 and 11 Ma (Sawyer et al., 1994) from large caldera complexes only~10 km to the north. Simulations indicate that these Miocene magma chambers would have disturbed local heat-flow regimes on the multi-million-year time scales producing elevated UZ temperatures to 6 Ma or younger Whelan, 2000, 2001;Whelan et al., 2001). (Whelan et al., 2002, p. 746-747) Modeling of thermal history indicates that the prolonged cooling of the UZ is consistent with heat flow produced by the regional magmatic activity responsible for widespread silicic volcanism 15 to 11 Ma, including the tuffs at Yucca Mountain. (Whelan et al., 2003, Sect. IV) Simulations of temperatures in the upper crust around a large intrusive heat source to the north indicated that cooling of the UZ at Yucca Mountain could have taken until 4-6 Ma, which is consistent with the fluid inclusion thermochronology . . . (Whelan et al., , p. 1884 These highest temperatures of deposition are explained in Marshall and Whelan (2001) with a thermal model. . . that links the slow cooling of the UZ at Yucca Mountain to the cooling magma body beneath the Timber Mountain caldera complex. . . . The driving mechanism is dissipation of a large amount of thermal energy emplaced at shallow crustal levels via magmatic processes. (Marshall et al., 2005, p. 221) Simulations using the HEAT code demonstrate that prolonged cooling of the unsaturated zone is consistent with magmatic heat inputs and deep-seated (sub-water table) hydrothermal activity generated by the large magma body~8 km to the north that produced the 15-11 Ma ash flows and ash falls that www.geosci-model-dev.net/7/1583/2014/ Y. V. Dublyansky: Evaluation of the US DOE's conceptual model of hydrothermal activity make up Yucca Mountain. (Whelan et al., 2008(Whelan et al., , p. 1041 Slow cooling of the unsaturated zone is, however, consistent with thermal models of regional heat flow around the large magma chamber~8 km north of Yucca Mountain that produced the voluminous 15-to 11-Ma ash flow and ash-fall tuffs that compose Yucca Mountain. The simulations approximate the thermochronology of secondary mineral deposition recorded in the unsaturated zone. (Whelan et al., 2008(Whelan et al., , p. 1072 Two-dimensional conductive/convective thermal modeling of the upper crust at Yucca Mountain indicates that modern geothermal gradients may not have been established at the repository horizon until ca. 6 Ma . These simulations are consistent with the cooling history derived from secondary minerals. Warmer saturated-zone groundwater, higher watertable elevations, and greater overburden thickness could have contributed to the prolonged period of unsaturated-zone cooling (Fig. 5 in Whelan et al., 2001). (Paces and Whelan, 2012, p. 245).

C2 Acceptance by DOE
Despite the lack of documentation and validation of the MICH thermal model, discussed in the previous section, the model was accepted by the US Department of Energy in 2001 and implicitly included in the Yucca Mountain Science and Engineering Report (DOE, 2001). In the aforementioned document, the elevated temperatures measured in fluid inclusions were dismissed by relating them to "a welldocumented thermal period that affected the volcanic rock for a long time after its formation" (p. 4-402) When technical data on the USGS modeling were published by the DOE contractor (BSC, 2004a) the lack of validation of the MICH model became apparent. Surprisingly, despite the numerous problems discussed in it, the document offered the following overall conclusion regarding the MICH model: "in summary, while the thermal model simulations of Whelan (2001, 2004) do not predict a thermal event that is as prolonged and pronounced as that recorded by secondary minerals at Yucca Mountain, their general model provides a mechanism to account for the presence of elevated temperatures between 10 and 6 Ma." (p. H-12).
The previous statement does not appear to be a sound scientific judgment. It is difficult to accept the notion that a model which failed, by a large margin, to match the benchmark empirical data "provides a mechanism" to account for these data.
In an analysis model report discussing features, events, and processes to be considered in the Yucca Mountain Total System Performance Assessment for License Application, the MICH thermal model was used as a key argument, to base the exclusion of hydrothermal activity from consideration (BSC, 2004b) (see analysis in Dublyansky, 2007). The latest (pre-License Application) version of this report (SNL, 2008) also relies on the purported positive outcomes of thermal modeling.
In 2004, DOE published a major document, Yucca Mountain Site Description (BSC, 2004c), which supports the Yucca Mountain License Application. This document also contained an inaccurate statement with respect to the validation of the MICH model: Two-dimensional conductive/convective thermal modeling of an idealized upper-crustal section at Yucca Mountain indicates that modern geothermal gradients were established by 6 Ma (Marshall and Whelan 2000). These results are consistent with the cooling history observed in unsaturated zone secondary minerals. (BSC 2004c, p. 7-81) Although Yucca Mountain thermal evolution based on these data require that heat dissipation from Miocene magmatic activity extended through longer periods of time than those expected by some , results are consistent with heat-flow models involving multiple injections of large magma bodies at shallow crustal levels. (BSC 2004c, p. 7-84)

C3 Acceptance by NRC
The US Nuclear Regulatory Commission staff reviewed the modeling results presented in the (BSC, 2004a) document. The review states: "additional thermal modeling by Marshall and Whelan (2001) suggests that the long-lived, nearsurface thermal perturbation at Yucca Mountain could not be reproduced by their thermal models, which predicted much faster cooling than inferred from oxygen and strontium isotope analyses in secondary minerals. . . " (NRC, 2005, p. 15).
The NRC review cites the four auxiliary hypothetical mechanisms proposed in BSC (2004a) to explain greaterthan-modeled temperatures and cooling times at Yucca Mountain (i.e., sustained magmatism, magmatic intrusions outside the caldera, additional overburden, and lateral outflow of hydrothermal fluids). The NRC report was explicit in stating that none of these hypothetical mechanisms is supported by site-specific factual evidence: There is little evidence to support DOE's scenario of sustained magmatism and, thus, sustained heating of crustal rocks within the Timber Mountain caldera. . . DOE does not cite any information to support the scenario of a hidden magmatic intrusion occurring south of the Timber Mountain or Claim Canyon caldera boundaries. DOE does not provide a technical basis to account for the thickness of potentially missing deposits needed for this scenario. In addition, DOE does not discuss how much additional burial would be needed in this scenario to account for paleotemperatures measured in 6 to 11 million year minerals at Yucca Mountain.
The NRC review rejects the first three mechanisms, but does accept the fourth: "subsurface outflow of hydrothermal fluids from the Timber Mountain caldera system, however, appears a credible scenario to account for elevated paleotemperatures preserved in 6 to 11 million year minerals at Yucca Mountain." (NRC 2005, p. 17).
The only argument put forth by the NRC staff to justify this acceptance was a statement that: ". . . such flows are commonly observed in geothermal systems that occur above and adjacent to large-volume magma bodies (e.g., Goff et al., 1988)" (p. 17). The reason why the NRC reviewers decided to cite general observations on geothermal systems and not consider the substantial body of site-specific information is unclear. Subsurface outflow of thermal fluids is known to have occurred at Yucca Mountain between 10.5-11.0 Ma, but only affected the deep-seated part of the rock sequence (below ca. 1000 m). Mineralogical, fluid inclusion, and isotopic evidence (Bish and Aronson, 1993, Sect. "Paleogeothermal conditions"; Feng et al., 1999, Sect. 4.1) and thermal modeling (Dublyansky and Polyansky, 2007, Sect. 4.3.3) indicate that this hydrothermal system did not cause heating of the unsaturated zone rocks commensurate with the fluid inclusion temperatures. If the site-specific information were used, this hypothetical mechanism would also have to be rejected. Acceptance of the "subsurface outflow" conjecture allowed NRC to accept the MICH model as a whole: Although studies of secondary minerals at Yucca Mountain by several organizations continue to this date, the NRC staff consider the conceptual model proposed by DOE for secondary mineral deposition at Yucca Mountain is generally consistent with available lines of evidence, notwithstanding remaining uncertainties in the age, timing, and origin of the thermal perturbations that produced elevated temperatures evidenced by fluid inclusions. . . (NRC, 2005, p. 17).
In view of the discussion above, this NRC decision to accept the MICH model does not seem to be scientifically justified.

C4 Acceptance by scientific community
The model of Marshall and Whelan (2001) was used to explain the elevated temperatures measured in the Yucca Mountain fluid inclusions by Wilson et al. (2003): "elevated temperatures within the sequence could be expected for a few million years following intrusion of the Timber Mountain Caldera at around 10 Ma (Marshall, 2000)." (Sect. 7.6). Bryan et al. (2009) used the MICH model to constrain the thermal history of Yucca Mountain in their study of the dissolution rates of feldspar in the Topopah Spring Tuff: "the available thermochronologic data for the Topopah Spring Tuff were compiled by Whelan et al. (2008). They also generated a thermal model, involving magmatic heat input from the nearby Timber Mountain volcanic center, which fit the thermochronologic data . . . " (Sect. 3).

C5 Summary
Having been cited in many technical publications over the past 12 years, the notion that the meteoric-infiltrationconductive-heating (MICH) model has been validated by thermal modeling has become firmly established in the scientific record. Meanwhile, examination of the available modeling results reveals no modeling outcomes which produce a reasonable match with the empirical benchmark data.