**2. Energy Balance**

To form a large magma body, the rate of heat input via magma injection must exceed the rate of heat loss via conduction, convection, and advection. For a one-dimensional system (e.g., a horizontal sheet of magma su fficiently broad that it can be considered one-dimensional; Figure 1), this is simply thermal energy in the form of magma put into the bottom of the system balanced against loss of heat out the top. The latter quantity must average out to surface heat flow, because heat loss is ultimately to the surface.

**Figure 1.** Schematic view of one-dimensional thermal model of a magma sheet with thickness much less than horizontal dimensions. Whether or not the magma body grows depends on the balance between heat in the bottom (as magma) and heat out the top. As heat eventually escapes to the surface, over the long term heat out the top of the sheet averages out to surface heat flow.

#### *2.1. Conductive Heat Flow*

Heat flow by conduction obeys Fourier's Law:

$$q = k \frac{dT}{dz} \tag{1}$$

where *q* is heat flux, *k* is thermal conductivity, *T* is temperature, and *z* is depth (positive downward). Surface heat flow in non-magmatic areas of the world is typically ~35–70 mW/m<sup>2</sup> [5], and the heat flux into the base of the crust is typically about half of this [6]. In this paper I use 40 mW/m<sup>2</sup> as a typical surface heat-flow value in non-magmatic continental areas. Given a typical *k* of 2 W/mK[7], a corresponding geothermal gradient, assuming conduction alone, would be 20 ◦C/km.

#### *2.2. Convection and the Nusselt Number*

Convection of fluids around a magma body can extract heat much more e fficiently than conduction alone. In geothermal areas convection is accomplished predominantly by circulation of meteoric water through fractures in rocks above and within a body of cooling magma. The ratio of total heat transfer by convection, advection and conduction to that of conduction alone is the Nusselt number, *Nu*:

$$Nu = \frac{\text{total heat flux}}{\text{convective heat flux}} \tag{2}$$

In geothermal areas, *Nu* can vary from unity for purely conductive systems to 100 or higher for systems with vigorous convection [8–10].

In magmatic areas regional heat flow can exceed 150 mW/m2. Individual thermal springs and fumaroles in the U.S. portion of the Cascade Range release up to several tens of MW over small areas [2] and collectively release ~1 GW, mostly in the Oregon and northern California segments. The Yellowstone hydrothermal system produces ~5 GW, corresponding to heat flows of ~2 <sup>W</sup>/m<sup>2</sup> averaged over the caldera system [11]. The far smaller (~100 km2) Grímsvötn geothermal area under Vatnajökull in Iceland produces comparable power, leading to the astounding heat flow estimate of ~50 <sup>W</sup>/m<sup>2</sup> [12].

Such high heat flows and power outputs indicate that the majority of heat transport is by convection of geothermal fluids rather than by conduction. Fournier [13] noted that the conductive geothermal gradient needed to sustain heat flow of 2 <sup>W</sup>/m<sup>2</sup> at Yellowstone would be ~1000 ◦C/km, requiring molten rock <1 km below the surface. This is contradicted by drill holes showing temperatures ~300 ◦C at 1 km [14] and by lack of seismic evidence for significant molten rock at such shallow depths [15]. Thus, high heat flow in hydrothermal systems such as Yellowstone are obvious evidence of the role that convection plays in moving heat.

## *2.3. Magma Accretion*

For a sheet of partially molten rock such as that in Figure 1, the energy balance between heat input into the bottom and released out the top can be approximated by a dimensionless magma accretion number

$$M = \frac{\rho LV}{\left(q - q\_{\text{big}}\right)}\tag{3}$$

where *q* is the heat flux out the top of the body, *qbkg* is the background heat flux, ρ is magma density, *L* is latent heat of crystallization, and *V* is the rate at which new magma is accreted to the body. For *M* < 1 new magma freezes as it arrives; for *M* > 1 the body of partially molten rock grows.

Using ρ = 2700 kg/m3, *L* = 3 × 10<sup>5</sup> J/kg, and *qbkg* = 40 mW/m2, we can estimate magma accretion rates that would be necessary for *M* > 1, and thus high enough to enlarge a growing magma body, as a function of surface heat flow, assuming steady-state (Figure 2). Although the relationship is linear, it is plotted on logarithmic scales to better show the relationship at lower *q*. Required vertical magma accretion rates are on the order of 1–10 mm/a (= 1–10 km/Ma) for some of the geothermal areas, but on the order of 100–1000 mm/a for the larger ones. Millimeter per year displacements of the Earth's surface are within the range of observed rock uplift and erosion rates (e.g., [16,17]), whereas the higher rates, involving injection of a pile of magma of thickness comparable to that of the continental crust in 1 Ma, are not. Geothermal areas typically have lifetimes on the order of 105–106 years [18], and sustained intrusion rates on the order of 100–1000 km/Ma are clearly incompatible with the geologic record.

**Figure 2.** Plot of values of accretion rate and surface heat flow for which the magma accretion number *M* = 1 (gray boundary). For accretion rates below this, cooling is fast enough that an incrementally emplaced pluton grows; for higher accretion rates, a body of partially molten rock grows. Black dots give estimated surface heat flow at selected geothermal areas, showing the magma accretion rates that would be needed to sustain those heat flows by conduction alone. The fields with high surface heat flow require magma accretion rates far in excess of geodetically observed rates.

#### **3. Cooling E**ff**ects of Circulating Fluids: Nusselt Number**

The Nusselt number *Nu* has a strong effect on *M* because it controls surface heat flow *q*. We can rewrite Equation (3) as

$$M = \frac{\rho LV}{\left(q\_s \cdot \text{Nu} - q\_{bk\text{g}}\right)}\tag{4}$$

where *qs* is surface heat flow under purely conductive conditions. *M* is greatly decreased for *Nu* > 1; for *Nu* = 10, *M* is roughly 9 times smaller than for *Nu* = 1. Thus, a Nusselt number on the order of 10–100 has a drastic effect on cooling of an accumulating magma body.

## **4. Thermal Models**

Thermal modeling provides insight into the effects of conduction versus convection on development of bodies of partially molten rock. Models appropriate to intermediate-composition (e.g., dacitic) magma were run using the Matlab partial differential equation solver; details of thermal modeling are in Appendix A. Enhanced heat flow owing to convection was simulated by increasing thermal conductivity *k* by a factor of *Nu* in the rock above the magma injections (e.g., [19]).

Results of a series of thermal models, with magma injected in sheets 50 m thick at time-averaged vertical accretion rates *V* of 5–20 mm/a and with *Nu* = 1 and 10, are presented in Figure 3. Models ran until 10 km of magma had accumulated. This is an unreasonable amount of magma to inject over a time on the order of 1 Ma, but the models were run this long to see if the temperatures stabilized. I monitored the maximum temperature just prior to each new injection as a proxy for thermal maturation of the system. For injections at 5 km depth (Figure 3a) and *Nu* = 1, at an injection rate of 5 mm/a, the maximum temperature leveled out at <600 ◦C at 1 Ma; 10 mm/a brought the temperature to ~750 ◦C after 1 Ma; and 20 mm/a brought the temperature to 900 ◦C, forming a significant body of eruptible magma, in 500 ka. However, with convection turned on (*Nu* = 10), none of the models approached

magmatic temperatures; at 20 mm/a, the temperature leveled out at <400 ◦C. Models run with injections at 10 km depth (Figure 3b) reached temperatures 50–100 ◦C higher, but again none reached magmatic temperatures when convection was turned on.

**Figure 3.** Thermal models of emplacement of 50 m sills at depths of 5 km (**<sup>a</sup>**,**<sup>c</sup>**) and 10 km (**b**,**d**), at rates of 5, 10, and 20 mm/a and Nusselt numbers of 1 and 10. (**<sup>a</sup>**,**b**) Maximum temperatures attained just before next injection of magma as a function of depth. Dotted lines at 700 and 850 ◦C are model solidus and 50 wt% melt, respectively; the latter is the approximate threshold for a mobile magma body. Models were run until 10 vertical km of magma was intruded. At 5 km, only the highest accretion rate produces a mobile magma body; at 10 km, both 10 and 20 mm/a models reach the 850 ◦C threshold. Heat flow (**<sup>c</sup>**,**d**) is far higher for the convective models with *Nu* = 10 than for conductive models with *Nu* = 1, reaching 1 <sup>W</sup>/m<sup>2</sup> or more for the faster accretion rates.

Surface heat flow is dramatically enhanced when *Nu* = 10 (Figure 3c,d). For example, at an accretion rate of 10 mm/a, surface heat flow by conduction levels out at ~300 mW/m2, whereas at *Nu* = 10 it reaches >1000 mW/m2.

#### **5. Climatic E**ff**ects on Volcanism and Hydrothermal Systems**

The magma accretion number arguments and modeling above predict that growth of large magma bodies is enhanced when heat transfer to the surface via convection is inhibited—that is, when *Nu* is small. The main factors that can inhibit convection include (1) lack of circulating water and (2) decreased permeability owing to precipitation of minerals in fractures or to particularly tight rock.

Although it has long been known that volcanism affects climate, an ever-growing body of work shows that the reverse is also true: climate can affect volcanism and hydrothermal systems over both short (seasonal) and long (glacial-interglacial) time scales. For example, a number of studies have shown a weak but persistent link between intensity of volcanism and waxing and waning of glacial cycles (e.g., [20–24]). These correlations have been tied to crustal stresses induced by changes in the mass of ice on land and water in the oceans. Climate and the state of the hydrologic system affect the intensity of hydrothermal activity; in the Kenya rift valley, ages of deposits from inactive geothermal systems correlate with periods of high lake level in the rift [25], and deposition of geothermal travertine in Italy, Turkey, and New Mexico correlates with global and regional paleoclimate [26–28]. At Yellowstone, the pressure increase caused by thick glacial ice apparently increased subsurface temperatures by shifting the boiling-point curve [29]. At shorter time scales, geyser periodicity in Yellowstone is weakly tied to both seasonal and decadal precipitation changes [30], and the power output of the El Tatio hydrothermal system in Chile, estimated from chloride flux, is twice as high in wet periods as in dry ones [31].

In many regions, fossil deposits indicate the presence of more vigorous and widespread geothermal activity in the geologically recent past. At the Coso field in California, extensive siliceous sinter deposits [32] underlie a 234 ka basalt flow that is unaffected by thermal activity, indicating more widespread activity in the Pleistocene [1]. Such extensive deposits of siliceous sinter are common in geothermal areas [33–36], and are typically interpreted as evidence of more widespread activity in the past, although migration of thermal vents over time can spread sinter over a wider area than the currently active one. Extensive thermal travertine plateaus and ridges [26,28] may also be relics of more widespread thermal activity in the past, presumably during wetter periods.

#### **6. Climate and the Growth of Persistent Magma Bodies**

Sufficiently rapid magma influx can always produce a large body of partially molten rock if the magma is trapped, but for systems near the critical balance between magma input and heat loss (Fig. 2), climate can play a role in whether the body attains large size because the rate of cooling depends on whether heat transport occurs by conduction or convection. In particular, in areas where precipitation is not sufficient to replenish water vapor emitted and lost during hydrothermal activity, convection can be starved, changing the mechanism of heat release from dominantly convective to dominantly conductive and dropping *Nu* from values of 10 or more to near unity. Thermal modeling (Figure 3) demonstrates that such a change can have a profound effect on the thermal state of an accumulating magma body, and that long-term dry climatic conditions should favor the development of caldera-forming magma bodies.

Coso exemplifies a geothermal field in an area that has undergone radical climate swings in the Quaternary. The field is currently quite arid, receiving ~200 mm/a annual precipitation, and is bordered on the south by China Lake playa, a desiccated remnant of a pluvial lake, on the north by Owens Lake playa, a former saline lake, and on the west by a dry valley once occupied by the pluvial Owens River. During the last glacial maximum, however, Owens Lake overflowed into the pluvial Owens River and fed a chain of lakes, including China Lake, that eventually spilled into Death Valley. Glaciers came to the valley floor in Owens Valley, which was occupied by tree species that indicate a far wetter climate [37]:

*It is suggested that meltwater from the retreating glacial ice inundated the Owens River Lake chain causing pluvial Owens Lake to reach its highstand. This caused an increase in e*ff*ective moisture, due to high groundwater, allowing the mesophytic Rocky Mountain juniper to exist at the site.*

#### **7. Implications of Precipitation Control on Magmatism**

#### *7.1. Paleogeographic Settings of Silicic Caldera Complexes*

The hypothesis that climate can affect magmatic systems by modifying the rate at which they cool may explain the paleogeographic settings of a number of large explosive silicic centers, as several of the world's largest silicic caldera provinces (Figure 4) formed in arid to hyperarid environments. The most obvious contemporary example of this is the Altiplano-Puna area [38] of the high Andes. Along the Andes major calderas are largely concentrated in the arid Central Volcanic Zone, with the greatest concentration in the driest part near the Atacama Desert (Figure 5). Extensive silicic ignimbrites of the

Early Cretaceous Paraná-Etendeka province of South America and Africa are intercalated with thick aeolian sandstones that indicate arid conditions throughout the eruptive sequence [39,40], consistent with the hot, arid, equatorial setting of central Gondwana during the Early Cretaceous [41]. Extensive Oligocene ignimbrites in Colorado, Utah, New Mexico, and Arizona are similarly intercalated with thick eolianites that are interpreted as remnants of a widespread erg [42] and corresponding arid conditions.

**Figure 4.** Locations of several large silicic ignimbrite provinces that formed in arid or plausibly arid environments, along with the three silicic large igneous provinces identified by Bryan [43]. All but the central Andean province formed during continental rifting, but most of these also formed in arid or hyperarid environments. Although there is a strong global correlation between magmatism and lithospheric extension, in detail extension may actually work to suppress growth of large magma bodies and caldera formation by opening fractures and promoting hydrothermal circulation.

Several other major silicic ignimbrite provinces formed in regions that were transitional from temperate to arid belts on paleoclimatic maps of Boucot et al. [44]. These include the Oligocene Sierra Madre Occidental province of Mexico [45–47] and explosive silicic volcanism of the contemporaneous Afro-Arabian province in Ethiopia and Yemen [48]. The latter accompanied eruption of the Ethiopian Traps during rifting that formed the Red Sea and Gulf of Aden. The Chon Aike/West Antarctica silicic province formed in the Jurassic during rifting of Gondwana [49,50]. The paleogeographic settings of these regions during explosive volcanism are not fully understood, but they were plausibly arid.

**Figure 5.** Locations of major late Cenozoic calderas (generally >10 km in diameter) in the Andes, from [51,52]. Map shows mean annual rainfall in mm from 190–2000, from worldclim.org. Calderas are dominantly found in the arid to hyperarid parts of the range.

#### *7.2. Role of Extension and An Explanation for Really Obvious Counterexamples*

The list of significant silicic caldera systems in non-arid environments is large and includes Toba [53], Krakatau [54], and numerous calderas in Japan [55] and New Zealand [56]. These clearly formed in regions with abundant water in the form of precipitation or seawater, but this observation does not negate the basic physics presented here. There are two obvious ways to promote development of an eruptible magma body in an area with abundant hydrothermal fluid. The first is to supply magma to the bottom of the system fast enough that even a vigorous hydrothermal system with high *Nu* cannot remove it fast enough to prevent solidification. However, modeling presented above indicates that extremely high magma accretion rates, on the order of 150 mm/a (= 150 km/Ma) are required to build up an eruptible magma body. This may be possible, but evidence for sustained crustal inflation at such extreme rates has not been recognized. This is a frontier area that deserves further investigation [57,58].

A second, more realistic way to increase *M* in a hydrothermal system is to seal up permeability so that convection cannot take place even when abundant working fluid is present. Self-sealing of hydrothermal systems is widespread and leaves behind clear physical evidence in the form of veins choked with precipitated minerals [59,60]. Wisian and Blackwell [61] noted that

*At least one major reason to associate geothermal systems with young* [extensional] *faulting is self-sealing, the process whereby cooling, ascending fluids precipitate minerals in pores and fractures (thereby reducing permeability). This process will eventually limit or eliminate flow in high-temperature cases. There is abundant evidence of self-sealing in many fossil geothermal systems* ...

Worldwide there is an undeniable, strong link between magmatism and lithospheric extension, but in detail there are anomalies in this relationship that sugges<sup>t</sup> a role for extension in suppressing development of large silicic magma bodies. By reasoning above, opening of sealed fractures is favored

by extensional faulting; thus, a system in which self-sealed fractures are episodically broken open by extensional faulting should have lower *M* than one in which fractures are simply sealed up, and systems without extensional faulting should be more prone to sealing, high *M*, and development of crustal magma bodies and large silicic eruptions. Analysis of the stress state of calderas around the world is beyond the scope of this paper, but, based on surface geology and earthquake focal mechanisms, both Sumatra and Honshu are undergoing shortening [62–64]. The North Island of New Zealand, which is rifting [65], is a major counterexample to the hypothesis proposed here.

The process of seal-breaking extensional faulting may explain the curious anticorrelation between extension and caldera formation in the Great Basin of the western United States, where eruption of major ignimbrite sheets tended to occur in between periods of extension and stratal tilting [66–68]. Gans and Bohrson [69] explained this by proposing that rapid extension can suppress volcanism by a number of mechanisms, including allowing access of meteoric water to midcrustal depths. Extension may suppress caldera formation by continually reopening permeability, promoting hydrothermal convection and magma body cooling. This process may also explain the paucity of calderas in arcs with backarc spreading [52].

#### *7.3. Possible Orographic Control of Caldera Formation*

Voluminous explosive silicic eruptions blanketed Nevada and western Utah in the Oligocene and Miocene during the southward-sweeping ignimbrite flareup [70,71], leaving behind 23 recognized calderas in Nevada and a nearly equal number of ignimbrite sheets that were probably erupted from unrecognized calderas (Figure 6). Magmatism swept southward across Nevada from ~45 Ma to ~15 Ma, but the caldera belt is bounded by a line that runs NNW, parallel to the Sierra Nevada Range, 100–200 km east of its crest. West of this line only the only recognized calderas are the 760 ka Long Valley caldera [72], 9.5 Ma Little Walker caldera [73,74], and 16 Ma Woods Mountains caldera [75].

It is possible that this sharp southwestern boundary to the caldera belt was caused by the westward transition from the arid interior West to the wetter continental margin in California. North America has been at roughly the same paleolatitude for the past 100 Ma [76], with the Great Basin in the belt of westerlies. Today, orographic precipitation wrings east-moving storms of their moisture over the coastal mountains and Sierra Nevada, producing the deserts of southern California, Nevada, and Arizona, including hyperarid Death Valley. In the mid-Cenozoic Nevada was likely a high plateau; regardless of whether the current Sierra Nevada rose above this plateau or were its western slope, the plateau was likely arid [77], promoting magma accumulation in the crust and caldera-forming eruptions.

**Figure 6.** Map of Nevada and surroundings, showing calderas from Henry and John [70] and mean annual precipitation. The ignimbrite flareup started at ~45 Ma in Idaho and swept south and southwest to southern Nevada ~15 Ma, but the voluminous caldera-forming events stopped as they approached the California border, perhaps fearful of high taxes or that their fruit would be confiscated. One explanation for this western border is that caldera-forming eruptions were favored in the arid high plateau of Nevada but suppressed in wetter regions closer to the ocean. LW, LV, and WM refer to Little Walker, Long Valley, and Woods Mountains calderas.

#### *7.4. A Final Speculation*

Cather et al. [77] proposed that large silicic eruptions can force world climate into a colder phase by fertilizing the oceans with iron, leading to carbon dioxide drawdown. Such a mechanism could interact with the processes discussed here in interesting ways. For example, if a warm climate leads to drier conditions in magmatic areas, volcanism could shift from effusive to explosive, triggering cooling that could dampen the climate and suppress such effusive volcanism.
