Research ArticleGEOLOGY

Geomorphic expression of rapid Holocene silicic magma reservoir growth beneath Laguna del Maule, Chile

See allHide authors and affiliations

Science Advances  27 Jun 2018:
Vol. 4, no. 6, eaat1513
DOI: 10.1126/sciadv.aat1513


Large rhyolitic volcanoes pose a hazard, yet the processes and signals foretelling an eruption are obscure. Satellite geodesy has revealed surface inflation signaling unrest within magma reservoirs underlying a few rhyolitic volcanoes. Although seismic, electrical, and potential field methods may illuminate the current configuration and state of these reservoirs, they cannot fully address the processes by which they grow and evolve on geologic time scales. We combine measurement of a deformed paleoshore surface, isotopic dating of volcanism and surface exposure, and modeling to determine the rate of growth of a rhyolite-producing magma reservoir. The numerical approach builds on a magma intrusion model developed to explain the current, decade-long, surface inflation at >20 cm/year. Assuming that the observed 62-m uplift reflects several non-eruptive intrusions of magma, each similar to the unrest over the past decade, we find that ~13 km3 of magma recharged the reservoir at a depth of ~7 km during the Holocene, accompanied by the eruption of ~9 km3 of rhyolite. The long-term rate of magma input is consistent with reservoir freezing and pluton formation. Yet, the unique set of observations considered here implies that large reservoirs can be incubated and grow at shallow depth via episodic high-flux magma injections. These replenishment episodes likely drive rapid inflation, destabilize cooling systems, propel rhyolitic eruptions, and thus should be carefully monitored.


The gradual accumulation of rhyolitic magma in the upper crust can promote large caldera-forming eruptions (1, 2), but the processes by which this occurs remain poorly understood (3, 4). A key issue is whether the magma flux into the upper crust is sufficiently large, over long enough periods of time, to sustain growth of magma reservoirs thermally capable of producing large rhyolitic eruptions, rather than crystallizing into plutons (58). Many of the observations fueling this debate come from radioisotopic dating of minerals [for example, (8)] or trace element diffusion clocks preserved within them [for example, (4)]. On the one hand, magma fluxes typical of plutonic systems are thought to normally be too low to sustain large, eruptible magma reservoirs in the upper crust (6, 912). On the other hand, a growing body of evidence suggests that silicic magma is stored long-term at relatively cool, nearly subsolidus conditions and is episodically remobilized by rapid injection of hot recharge magma that may propel destabilization and eruption (3, 4, 8, 1315).

The magnitude and pattern of surface deformation offer another important means to understand magmatic processes operating beneath restless, occasionally active, volcanoes (16, 17). Measurement of deformation by satellite geodesy in caldera volcanoes that produced supereruption scale (18, 19), as well as modest volume (2023) rhyolitic eruptions, has revealed inflation affecting regions of hundreds of square kilometers, over periods of months to years, typically at rates of 10 cm/year or less (17). Arrival of new magma into an extant, shallow magma reservoir is the common explanation for the surface inflation [for example, (22)], although the pressurization of magmatic fluids also contributes in some systems (17, 24). What has remained elusive is a means of measuring a long-term flux rate that integrates several magma pulses and that may leverage interpretations of whether a magma reservoir is likely to grow and erupt or to freeze into a pluton (7).

The Laguna del Maule (LdM) volcanic field (Fig. 1) comprises the greatest concentration of postglacial (younger than ~20,000 years) rhyolite in the Andes and includes the products of ~40 km3 of explosive and effusive eruptions (2528). Recent observations at LdM by interferometric synthetic aperture radar (InSAR) and global positioning system (GPS) satellite geodesy have revealed inflation at rates exceeding 20 cm/year since 2007 (2931), thereby capturing an ongoing period of growth of a potentially large upper crustal magma reservoir (27). The current episode of inflation has been explained by a model of transient supply of magma into this reservoir at a depth of 4.5 km and requires recharge at a rate of 0.03 to 0.04 km3/year (31, 32), which is a flux sufficient to destabilize a cool silicic magma reservoir (6, 7, 14, 22). Here, we use a geomorphic record of surface deformation at LdM that offers an unprecedented opportunity to link this current episode of unrest and inflation in a rhyolite-producing system to the record of rhyolitic volcanism and magma intrusion spanning the last 10,000 years.

Fig. 1 Maps of postglacial lava flows in the LdM volcanic field [modified from (25, 34)].

Shown are postglacial rhyolites (pink), rhyodacites (orange), and andesites (pale green). The presumed vent for the plinian eruption of the Rhyolite of LdM (unit rdm) is beneath the lake. The star in the inset shows the location of LdM in South America. Important map units include (i) rle, Rhyolite of Los Espejos that dammed the lake at 19 thousand years; and (ii) rcb, the many rhyolitic lavas of the Barrancas complex that were emplaced during the last 14 thousand years. These and other map unit abbreviations follow those in (25, 34).


Rhyolitic volcanism and unrest

LdM sits atop the southern Andean range crest where, following a rapid retreat of glaciers between ~23 and 19 thousand years ago (based on 40Ar/39Ar dating of pre- and postglacial lava flows), a flare-up of dominantly silicic volcanism occurred within the central lake basin (Fig. 1) (25, 33, 34). The silicic eruptions include effusive and explosive events and are volumetrically dominated by crystal-poor rhyolite. Rhyodacite and andesite eruptions also occurred throughout postglacial time but are concentrated in the northwestern to western LdM basin several kilometers from the locus of rhyolitic volcanism (Fig. 1). Rhyodacitic lavas are distinguished from rhyolitic ones by higher crystallinities, common amphibole crystals, and abundant centimeter-scale chilled mafic inclusions, which are nearly absent within the rhyolites (25, 33, 34). 40Ar/39Ar dating and volcanic morphology indicate that the frequency and spatial extent of the rhyolitic eruptions were greatest during two pulses: an early postglacial period (22 to 19 thousand years ago) and the middle to late Holocene (Fig. 1) (33, 34). The early postglacial rhyolite eruptions began with the ~20-km3 plinian Rhyolite of LdM (unit rdm) (28), the large volume of which distinguishes it from the subsequent, smaller-volume (0.5 to 3 km3) rhyolites. Between 14.5 and 8.4 thousand years ago, rhyolite erupted only from the Barrancas complex in the southeastern lake basin, after which the focus of rhyolite volcanism expanded northwestward and northward (Fig. 1) (25, 33, 34), while also continuing within the Barrancas complex (35). Excluding ash fall deposits (28), the minimum total volume of rhyolitic lava and pyroclastic flows that erupted from at least 10 vents ringing the southeastern margin of the basin during the last 14.5 thousand years is ~8.4 km3, more than half of which erupted during the last 5.6 thousand years (33). This chronology is supported by 36Cl surface exposure ages (Materials and Methods) of 2.5 ± 0.7, 1.4 ± 0.6, and 0.8 ± 0.6 thousand years from three late Holocene rhyolite flows that have little ash cover and are thus among the youngest eruptive features in the basin (Fig. 2 and table S1A).

Fig. 2 Map of preserved highstand paleoshoreline (color-contoured for elevation) and postglacial silicic lava flows.

Map unit abbreviations follow those in (25, 33, 34). ka, thousand years ago (ka ago).

Amphibole barometry, trace element compositions, thermodynamic modeling of phase equilbria, and Sr, Pb, and Th isotope ratios indicate that crystal-poor rhyolite erupted at LdM was produced in the upper crust via fractional crystallization and assimilation of previously emplaced magma ± juvenile crust at depths of 4 to 6 km (33, 34). Each cubic kilometer of erupted rhyolite is distilled from ~3 to 4 km3 of andesitic to rhyodacitic magma, suggesting that a substantial magma reservoir underlies the LdM basin (33). Subtle, temporally coherent variations in whole-rock trace element, Fe-Ti oxide, plagioclase, and zircon compositions indicate that the rhyolites evolved within at least two physically discrete reservoirs that may or may not have been contiguous beneath the LdM basin (33, 34). Mafic to intermediate eruptions are conspicuously absent within the locus of rhyolite volcanism. Yet, quenched basaltic andesite inclusions are common in the rhyodacitic lavas (Figs. 1 and 2), indicating that the ascent of mafic magma to shallow depths is frequent but is intercepted by a voluminous shallow silicic magma reservoir.

A Bouguer gravity survey was conducted in the LdM basin, with a focus over the recent area of inflation (36). Spatially coincident with the center of current uplift is a negative Bouguer gravity anomaly that has been modeled as a low-density, rhyolite-bearing magma reservoir (36). Inversion of the density distribution yields a model wherein a 30-km3 body of low-density rhyolitic magma containing exsolved fluids at depths of 2 to 4 km is embedded within a 115-km3 domain comprising denser crystal mush with several-percent melt (36). Rhyolite is extracted from rhyodacite magma near the base of the low-density domains inferred from the gravity field (33, 34, 36). Underlying this shallow reservoir is a long-lived transcrustal magma system [for example, (3, 37)] that likely extends downward into a deep crustal zone of magma storage, mixing, and assimilation (33, 38).


Dating the highstand paleoshoreline

The Espejos rhyolite coulée, 40Ar/39Ar-dated at 19.0 ± 0.7 thousand years before present (33) (Fig. 2, unit rle), dammed the northern outlet of LdM and raised the lake level by ~200 m to form a prominent basin-wide highstand shoreline (25). The lava dam was breached, resulting in an outbreak flood of 15 km3 of water pouring through the gorge and scouring the upper Maule river valley (25). The highstand shoreline is variably preserved and expressed around the basin. Atop the Espejos rhyolite coulée, it is a beach bar deposit of rounded boulders of the rhyolite (Fig. 3A). To the northeast, it is a bench up to several tens of meters wide cut directly into Pleistocene volcanic rocks (Fig. 3B). Around the eastern and southern edges of the basin, it is a bench cut into older Pleistocene lavas, covered by rounded beach boulders and late Holocene tephra (Fig. 3C). To determine the age of the highstand shoreline, we collected samples for surface exposure dating from five wave-cut terrace outcrops in the northern LdM basin at elevations between 2374 and 2388 meters above sea level (masl) (Fig. 2). Cosmogenic 36Cl measurements (Materials and Methods) yield ages of 9.4 ± 0.4, 8.8 ± 0.6, 7.5 ± 0.3, 6.6 ± 0.6, and 4.2 ± 0.2 ka (±2σ analytical uncertainties; table S1A). Two outcrops, one in rhyodacite ignimbrite at 2376 masl (Fig. 3B), the other 4 km to the northeast in andesite at 2386 masl, comprise tens to hundreds of square meters of block-jointed rock in windswept locations and yield the oldest dates of 9.4 ± 0.3 ka and 8.8 ± 0.6 ka that are indistinguishable from one another at the 95% confidence level. The three other sites on the rhyodacite ignimbrite spanning the same range in elevation yield younger ages. For samples that contain low native chlorine (table S1C), erosion of the surface or burial under ash, snow, or water will lower apparent 36Cl surface exposure ages, whereas for samples with high Cl, the impacts of erosion and burial tend to increase the apparent exposure ages. We therefore infer that the two younger ages from the shoreline derived from samples with low Cl (7.5 ± 0.3 thousand years ago and 4.2 ± 0.2 thousand years ago) reflect erosion or burial and that each gives a minimum duration of exposure. Likewise, the oldest 36Cl exposure age from the shoreline (9.4 ± 0.3 thousand years ago), also derived from a low-Cl sample, provides a minimum age of exposure. However, we consider this figure to be a closely limiting minimum age because posited rock surface erosion rates have a negligible impact (table S1A), and past cover by ash or snow is deemed minimal because of the windswept location of these outcrops. The exposure ages derived from moderate- to high-Cl material (8.8 ± 0.6 thousand years ago and 6.6 ± 0.6 thousand years ago) are more difficult to evaluate with regard to erosion, past cover, and pore water content (table S1A) but may in part reflect episodic submergence associated with seasonal changes in paleolake level similar to that observed along the shore of the modern lake (25). We conclude that the abrupt drop in lake level to near its modern elevation occurred at ca. 9.4 thousand years ago and left the highstand shoreline outcrops exposed to buildup of 36Cl. This conclusion is bolstered by the ages of two lavas that flowed onto the highstand shoreline. A surface on the Colada Dendriforme rhyodacite flow (unit rcd; Figs. 2 and 3D) yields a 36Cl age of 8.4 ± 1.2 thousand years, whereas a rhyolite flow on the north flank of the Barrancas complex (Fig. 3C) yields an 40Ar/39Ar age of 5.6 ± 1.1 thousand years (33).

Fig. 3 Outcrops illustrating geomorphic expression of highstand paleoshoreline in select localities around the LdM basin.

(A) Rounded boulders of 19.0 thousand years old Espejos rhyolite deposited on an angular outcrop of the same unit to define the highstand of the paleoshoreline (dashed yellow line). (B) Wave-cut terrace on 990 thousand years old Bobadilla rhyodacite ignimbrite north of LdM. The block below the hammer yields a 36Cl exposure age of 9.4 ± 0.4 thousand years. (C) Beach deposit of rounded boulders of Pleistocene andesite partially buried by rhyolitic ash and lapilli at the southeastern end of the lake basin (Fig. 2). The rhyolite lava from the Barrancas complex flowed over this shoreline terrace at 5.6 thousand years ago. (D) View north of the southern margin of the rhyodacite of Colada Dendriforme (unit rdcd in Figs. 1 and 2) where it flowed over a highstand paleoshoreline terrace (arrows). A glassy tower on the levee of the flow lobe on the right yields a 36Cl exposure age of 8.4 ± 1.2 thousand years that provides a lower limit for the age of the paleoshoreline. Photo credits: B.S.S. [taken on 1 January 2016 (A), 16 January 2015 (B), 8 January 2015 (C), and 29 March 2015 (D)].

Measuring Holocene surface deformation

We measured GPS positions of 64 sites on the highstand surface in a rapid static mode using five continuously recording GPS stations as base stations (Fig. 2), yielding elevations with centimeter-scale precision (Materials and Methods and table S2). We corroborated these GPS measurements by mapping each onto a photogrammetrically generated 1-m-resolution digital elevation model (DEM; Materials and Methods) and outlining the preserved highstand shoreline surface outcrops (Fig. 2 and fig. S1). Most surfaces are inclined toward the modern lake, and many are covered by tephra (Fig. 3C); hence, the main uncertainty in the GPS measurements comes from the choice of the GPS site location within one inclined “outcrop,” which may extend over 50 to 100 m in a horizontal distance. The comparison between the DEM and GPS height measurements allows us to estimate a total uncertainty for each GPS position of ±3 m (Materials and Methods and fig. S2). Elevations of the 64 GPS sites show 62 m of relief and range from 2366 to 2428 masl, with the lowest values along the western and northern perimeter of the lake basin and the highest along the south and east (Fig. 4 and fig. S1A).

Fig. 4 Model of single magma intrusion to explain Holocene surface deformation at LdM.

(A) Map of interpolated elevations of the paleoshoreline (color scale in meters) relative to the lowest measured level to the northwest (NW) (2372 m), as measured at 64 GPS sites (black circles) on the highstand paleoshoreline around LdM (blue outline). The white circle indicates the location of the best-fit prolate spheroid deformation source estimated from inversion of the GPS data set. (B) The analytical solution for a pressurized prolate spheroid source (54) has seven free model parameters: location (x0, y0), depth (z0), dimensionless pressure change from which we calculate the total volume change (ΔV), spheroid aspect ratio (A = b/a), dip angle (θ, positive downward), and strike angle (ϕ, clockwise from north). Gray histograms show the initial random distribution and bounds for each parameter. Blue histograms show the best-fit values for each parameter and each of the 250 runs. Final best-fit value (minimum χ2) for each parameter and uncertainty estimates from Jackknife resampling are indicated on top of each panel. (C) Vertical displacement data (red error bars) and modeled values (blue squares) as a function of the radial distance from the source.


The highstand paleoshoreline is inferred to have been a horizontal surface when the lake drained at 9.4 thousand years ago. Its elevation gradient and spatial distribution are inconsistent with predicted isostatic rebound in response to draining of the lake [for example, (39)] (Materials and Methods and fig. S3). The LdM basin is adjacent to the Andean range crest with current elevations of 2100 to 3000 masl throughout the region shown in Fig. 1. Here, during the Last Glacial Maximum ~21 thousand years ago, an ice cap ~80 km across and up to several hundred meters thick occupied the deepest valleys and range crest basins (40, 41). Deglaciation of the LdM basin was likely completed by 19 to 18 thousand year ago, resulting in isostatic adjustment to less than 5 MPa (megapascals) equivalent of lost ice. Modeling of glacial isostatic uplift of the European Alps (42), where the extent and thickness of ice at the Last Glacial Maximum were much larger than that at LdM, indicates that deglaciation induced about 2 mm/year of uplift. Applying this estimate as a maximum for LdM implies that less than 20 m of uplift, across the entire region shown in Fig. 1, may have occurred during the last 9.4 ka. This is inconsistent with both the short wavelength and the >60 m magnitude of surface deformation measured at LdM since 9.4 thousand years ago (Materials and Methods). Instead, the gradient is consistent with deformation associated with an inflating magma reservoir located below the southeastern portion of the lake basin near the Barrancas volcanic complex (Fig. 2). Moreover, the 60+ m of permanent deformation during the last 9.4 thousand years favors intrusion and solidification of magma, rather than pressurization of volatiles in a hydrothermal system.

To investigate the origin of the deformed paleoshoreline, we model the cumulative vertical displacement measured over the last 9.4 thousand years and relate it to the numerical modeling results of magma injection into a long-lived reservoir on the decadal time scale (32) The final GPS elevation data set is referenced in space to the lowest elevation in the northwest (Fig. 4) and inverted for three different deformation source geometries. Assuming that the crust is an isotropic, homogeneous, elastic half-space, we find the best-fitting source to be a prolate spheroid located in the southeast of the lake basin experiencing a volume change of 13 km3 over 9.4 thousand years at a depth of 7 km (Figs. 2 and 4A). A statistical F test confirms that this model with seven adjustable parameters is a significantly better fit than a simple spherical source with 95% confidence (Materials and Methods). The inflating source model reproduces the main characteristics of the spatial pattern and magnitude of the paleoshoreline deformation. The location of the inferred source of magma intrusion is centered among the rhyolitic volcanic complexes that ring the southern and eastern flanks of the LdM basin and from which eruptions were focused during the last 9.4 ka (Figs. 1 and 2).

Whereas there are no obvious offsets of the highstand paleoshoreline bench by faults, there is a relatively abrupt change in elevation of 20 m over a distance of 500 m across the inferred trace of the regional Troncoso fault (Fig. 2 and fig. S4). Cumulative movement of several meters along this fault may explain some of the larger residuals visible on the radial profile (Fig. 4C).

The kinematic model considers a constant change in pressure with time to reproduce the pattern of cumulative uplift. In the absence of a record of the temporal evolution of the uplift over the 9.4 thousand-years time interval, we use information revealed from the study of the ongoing uplift episode that began between 2004 and 2007. To reproduce the exponentially increasing and decreasing rates of uplift measured at LdM, Le Mével et al. (32) modeled the flow of magma through a conduit into an ellipsoidal source underlying most of the central part of the lake basin. If several such non-eruptive deformation episodes with the same characteristic time constant occurred over 9.4 thousand years ago, then 16 distinct episodes such as this would be necessary to uplift the paleoshoreline by 60 m as plotted in Fig. 5. However, the modeling results indicate two main differences between the deformation sources: The location has moved northwestward since 9.4 thousand years ago (Fig. 2), and the Holocene deformation source was deeper, leading to a wider deformation footprint at the surface.

Fig. 5 Sketch illustrating a conceptual model of episodic deformation at LdM.

Triangles on the upper x axis indicate the eruption ages of the rhyolite flows mapped in (25, 33) that issued from the Nieblas, Barrancas, Divisoria, and Cari Launa eruptive centers. Red filled triangles indicate flows dated by 40Ar/39Ar or 36Cl in (33) and this study, and gray filled triangles indicate ages inferred from field relationships; these flows have been placed at regular intervals between the bounding ages. The black square shows the cumulative vertical displacement estimate of 60 m over 9.4 thousand years. If the uplift rate had been constant over this time interval, then the long-term average rate would be 6.4 mm/year (dashed black line). The inset shows the temporal evolution of the vertical displacement as observed during the uplift episode that is ongoing since 2007 (red line) (32) and the projected rate assuming no further perturbations (dashed line). If all the deformation episodes share the same temporal and spatial characteristics, then 16 of these 50-year-long episodes would be needed to reproduce the 60-m total uplift of the paleoshoreline over the Holocene, with a magma intrusion average recurrence interval of 624 years.


Geomorphic constraints on magma reservoirs

Geomorphology has been used to constrain the depths and growth rates of magma reservoirs in a variety of settings, a few examples of which are outlined here. Whereas InSAR data show that growth of the magma body at a depth of 19 km below Socorro, New Mexico, propels recent uplift at 2.5 mm/year, geomorphic evidence suggests that it is unlikely that this uplift has persisted for more than a few centuries (43). Geodetically measured uplift rates at Uturuncu and Lazufre volcanoes in the Central Andes are between 1.0 and 3.5 cm/year (44). Several lines of geomorphic evidence suggest that at Lazufre, expansion of a magma reservoir at depths of 10 to 20 km generated pulses of surface uplift that total several hundred meters during the last 400 thousand years, punctuated by long intervals of quiescence or subsidence (45). Also in the Central Andes, a long-wavelength topographic dome about 1 km high reflects growth of the enormous 10- to 20-km-deep Altiplano-Puna magma body during perhaps the last 10 million years (Ma) (46). Resurgent doming or uplift of tens to hundreds of meters over periods of centuries to several millennia at rates of between 1 and 20 cm/year has also been ascribed to magma intrusion beneath calderas at Toba (47), Campi Flegrei (48), and Iwo Jima (49). Features that distinguish LdM from these other examples include the following: (i) the shallow depth of the magma reservoir inferred from gravity (36), geodesy (3032), and petrology (33, 34); (ii) a highly resolved series of rhyolitic eruptions of well-constrained volumes before, and during, the 60+ m of surface deformation of the last 9.4 thousand years (Fig. 1); and (iii) the exceptional magnitude of the current uplift at more than 20 cm/year during the last decade that is centered among the many rhyolitic vents (Fig. 2). We next integrate these observations to explore the long-term dynamics of the LdM magma reservoir.

Incremental growth and eruption of a shallow plutonic reservoir

Deformation of the paleoshoreline, coupled with knowledge of the volcanic output, the numerical simulations of magma intrusion, and geometry of the Bouguer gravity anomaly (33), leads to the model outlined in Fig. 6. During the last 9.4 thousand years, the ratio of rhyolite erupted (~9 km3) to the volume of magma transiting into and out of the reservoir (22 km3) is about 0.4, suggesting reservoir growth through underaccretion of magma as suggested in some epizonal plutons (6, 50, 51). Thermal modeling suggests that at a depth of ~7 km, the time-averaged magma flux here of 0.0023 km3/year is unlikely to favor the eruption of rhyolite, instead leading to solidification of a pluton (7, 10, 52). However, the >1.5 million years of LdM volcanism (25) attests to a magma flux that has thermally primed the upper crust to host a large, dynamic, eruptible reservoir of silicic magma [for example, (11, 33, 34)]. During the Holocene alone, at least 10 eruptions of crystal-poor rhyolite reflect crystallization of perhaps 40 km3 of andesitic to rhyodacitic magma within the shallow reservoir (34).

Fig. 6 Schematic evolution of the upper crustal magma reservoir beneath LdM during the last 22 thousand years.

The current Bouguer gravity low (36) is interpreted to include a smaller melt-rich body embedded within a larger, mostly crystalline, body centered under the lake. The gravity data suggest that the magma reservoir has begun to solidify into relatively dense crystal-rich mush beneath the southern end of the lake basin. Petrologic data (33, 34) imply that low-density rhyolitic melt segregates from crystal mush at depths of 4 to 6 km. The source of current inflation is best modeled as a sill opening at a depth of 4.5 km (31, 32). Thus, in the lowermost panel, we show that the partially molten zones and locus of magma injection coincide with the deeper portion of the inferred source of the Bouguer gravity low.

The absence of measured deformation before 2007 (29, 30) demonstrates that the current inflation event reflects an episodic process operating beneath LdM. Following the uniformity principle, the long-term deformation measured here would also reflect successive uplift episodes. We propose that the magma supply is best understood as an aggregation of repeated high-flux episodes exemplified by the ongoing inflation. During one such episode, the transient uplift rate could be as high as 200 mm/year, corresponding to a magma recharge rate of as high as 0.04 km3/year (Fig. 5) (32). These recharge episodes each propelled destabilization of the reservoir via the addition of volatiles and ascent of bubble-rich plumes of rhyolitic melt through low-porosity crystal mush (Fig. 6) (33, 53, 54). Moreover, the locus of magmatic intrusion, wavelength of surface inflation, and locations of rhyolitic eruptions have migrated several kilometers during the last 20 thousand years, highlighting the spatially heterogeneous, incremental growth of the magma reservoir. The Bouguer gravity anomaly implies that the reservoir today is partially molten under the lake (~1800 to 2400 kg/m3), with much denser crystal mush or subsolidus material (~2700 kg/m3) underlying the southern flank of the basin (36), suggesting that here it has cooled and begun to solidify (Fig. 6). Additional gravity measurements on the Barrancas complex and to its south are needed to test and refine this inference. The time-averaged rate of magma accumulation beneath LdM during the Holocene of 0.0023 km3/year is equivalent to many well-dated plutons (7, 9, 10, 50, 52) but has likely been punctuated by numerous high-flux recharge events at rates 20 times faster. Our findings illustrate that the incremental growth of plutons comprising several hundred cubic kilometers, accompanied by eruption of tens of cubic kilometers of rhyolite, is not restricted to deep hot crust. Rather, plutons may grow in the uppermost crust via repeated high-flux additions of hot magma to a gradually expanding, but mostly frozen, body of plutonic mush (4, 14, 34). Future explosive eruptions of modest to large volume are to be expected as a consequence of rapid recharge events akin to that observed today at LdM.


36Cl dating

For 36Cl exposure age determinations, five samples were collected from horizontal surfaces of the lake highstand shoreline notched into 990 thousand-years-old welded rhyodacite tuff or 1.01-million-years-old andesite north of LdM, and four samples came from horizontal surfaces on rhyodacitic or rhyolitic lava flows that exhibit uneroded morphologically youthful flow structures (Figs. 2 and 3 and table S1A). Local topographic maxima were selected for sampling to minimize the effects of shielding, which was evaluated using an inclinometer. Samples were collected using a hammer and chisel at least 5 cm from adjacent vertical surfaces to avoid edge effects in calculating ages. Whole-rock separates, crushed to 125 to 250 μm, were prepared from the uppermost 1.75 to 3.0 cm of each sample and ultrasonically cleaned in deionized water. Cl was extracted in clean laboratory facilities at the University of New Hampshire following methods developed by Stone et al. (55) and modified by Licciardi et al. (56). Crushed samples were pretreated with 2% HNO3 and spiked with an enriched 37Cl tracer and then dissolved in HF-HNO3 solution. Upon complete digestion, insoluble fluoride compounds were removed by centrifugation, and Cl was precipitated as AgCl with the addition of AgNO3. The precipitate was further purified by redissolution in NH4OH and the addition of BaNO3 to precipitate sulfate as BaSO4. AgCl was then reprecipitated by the addition of 2 M HNO3 and AgNO3, washed repeatedly in deionized water, and dried.

The 35Cl/37Cl and 36Cl/37Cl ratios were measured at the Center for Accelerator Mass Spectrometry at the Lawrence Livermore National Laboratory (LLNL) facility. Procedural blanks contributed 0.1 to 5.0% error to the 36Cl concentration errors in the unknowns, and appropriate blank corrections were made before age calculations. Ages were calculated using the online CRONUScalc 36Cl exposure age calculator and the LSDn scaling scheme (5759). Sensitivity analyses were conducted using the CRONUScalc calculator to isolate the potential impacts of rock surface erosion rates and pore water content on the apparent exposure ages (table S1A). A prescribed rock surface erosion rate of 5 mm per thousand years, which represents a plausible value determined on similar volcanic lithologies in the Quelccaya region of the Peruvian Andes (60), was seen to have a minimal effect on apparent exposure ages derived from samples with low native Cl [15-SLM-05, 22.9 parts per million (ppm) Cl; SLM-16-28, 25.3 ppm; SLM-16-29, 26.2 ppm], including the sample that provides the closely limiting 9.4 thousand years minimum age of the highstand shoreline. Likewise, a prescribed 0.5 fraction of pore water content had negligible impact on ages derived from these low-Cl rocks. In contrast, the same prescribed surface erosion rate and pore water content were seen to have variable and significant impacts on the apparent exposure ages of surfaces derived from rocks with moderate to high Cl contents (table S1, A and C).

Static GPS measurements

GPS data that measure shoreline height were collected during two survey periods in 2015 and 2016 using two Trimble dual-frequency receivers: 5700 and NetR9 with Zephyr and Zephyr Geodetic antennas. Each GPS point was occupied for at least 10 min, following a rapid static method. To obtain precise coordinates for each point, we followed the GPS differential procedure (61), using, as reference, five continuous GPS (cGPS) stations operated by the Observatorio Volcanologico de los Andes del Sur (OVDAS) located around the lake [Fig. 2; base stations are MAU2, PUEL, LDMP, COLO, and NIE2 cGPS sites of (31)]. Precise coordinates of continuous stations were estimated using Trimble postprocessing service (Topcon) relative to the World Geodetic System 84 (WGS84) ellipsoid. These positions were reprocessed with Trimble Business Center and compared with results from long time series processed with GIPSY-OASIS software (release 6.3) from the Jet Propulsion Laboratory, yielding similar results. For each campaign, 2015 and 2016, this procedure ensured that ongoing surface deformation is accounted for. Each GPS point was processed by calculating baselines with respect to the nearest or the two nearest permanent-fixed stations (fig. S2). We processed baselines and obtained fixed solutions for 64 GPS points; when two baselines were processed for the same point, we adjusted both results for a robust solution. Results are shown in table S2. Most GPS sites were chosen at the midpoint of the inclined slope of the paleoshoreline surface or at the highest break in its slope. The mean difference between the elevations in the DEM and the GPS sites was 0.27 ± 1.07 m. Of the 64 GPS sites, 48 (75%) agreed within 1 m with the DEM. The minimum and maximum differences were 0.02 and 2.62 m, respectively. Thus, we attributed a total uncertainty for each GPS position of ±3 m (fig. S2).

Digital elevation model

A DEM was constructed using photogrammetric analysis of several thousand oblique, near-vertical, digital images acquired during two helicopter survey campaigns in 2017 and methods outlined in (62, 63). Ground control points used in the modeling include several purpose-placed markings and several of the seismic instruments deployed to study the LdM magma system. The DEM has a spatial resolution of 1 m.

Isostatic adjustments to lake-level lowering and deglaciation

To consider surface displacement due to the sudden drop in lake level at 9.4 thousand years ago, we applied an elastic formulation (64) to calculate the unloading due to a 200 m drop in the water table that followed the current shoreline. Even an extremely low value of crustal shear modulus (1 GPa) would only account for 14 m of uplift, limited to the lake basin with a maximum centered within the lake (fig. S3). Therefore, we conclude that this process does not explain the observed pattern of uplifted shoreline to the SE by 60+ m (Fig. 4A).

To calculate the amount of uplift caused by viscoelastic postglacial rebound, we considered a model for the European Alps (42). Their estimate of 2 mm/year in the Alps was likely much faster than what the same process would produce anywhere at LdM in the Andes. First, the Alps sustained a larger (800 × 250 km) and much thicker (~2 km) ice cap. Second, since the Andes comprise a magmatically active subduction zone, the lithosphere is likely hotter and thus less viscous than below the Alps. The latter effect would reduce the effective elastic thickness and thus the contribution to uplift by deglaciation relative to mantle flow (42). In the Alps, the observed glacial uplift spans spatial length scales of hundreds of kilometers (42). At LdM, however, the warping of the LdM highstand paleoshoreline occurs over a much shorter spatial length scale of the order of ~15 km.

Both of these unloading processes (melting glaciers or draining lake) would lead to a rate of uplift that is fairly constant over the time scale of a decade. Neither process can explain the uplift measured between 2007 and 2017 at LdM. As discussed previously, the rate of uplift between 2003 and 2004 is 0 ± 10 mm/year (2931). We conclude that the doming of >60 m is instead driven by magma intrusion below the southern lake basin consistent with the eruptive history (Figs. 1 and 2).

Analytical modeling and interpretation

Before modeling, elevations were referred to a site in the NW assumed to be the original shoreline height. It is the furthest point from the source and has an elevation of 2372 m. The relative heights are the input for the analytical models. The horizontal displacement was unknown and was therefore set to 0 ± 30 m in the model input. All the models considered in this study assumed the crust to be an isotropic, homogeneous, elastic half-space. We tested three possible source geometries: (i) a sphere, implemented using the formulation of (65); (ii) a prolate spheroid (66); and (iii) a sill or penny-shaped crack using the three-dimensional Green function proposed by Fialko et al. (67). The inversion for the three different deformation sources was implemented in the dMODELS software package (68). The nonlinear inversion algorithm is a combination of a local optimization (interior-point method implemented in the MATLAB function fmincon) and a random search. For each set of parameters, 250 runs were calculated, and the parameter uncertainty estimates were then computed using a Jackknife resampling method.

We used analytical solutions of a simple deformation source in an elastic half-space to estimate the source volume change at depth from the cumulative surface uplift. Numerical models of the well-constrained temporal evolution of uplift for the ongoing episode of deformation started in 2007 helped us interpret the results in terms of magma volume estimates. Note that a magma injection model with a viscoelastic rheology for the crust yields up to 4% more displacement than the purely elastic model (32) and therefore would lead to a smaller estimate of magma injection volume for each episode. Considering the magma compressibility would, on the other hand, significantly affect the volume estimates of the magma injection. As shown in (32), the magma injection volume needed to explain the current uplift is twice as large if we assume a compressible magma (with B = 2.1 × 10−9 Pa−1). Consequently, the incompressible magma volume estimates, including the 12.7 ± 0.9 km3 prolate spheroid in Fig. 4B, are minimum estimates.


Supplementary material for this article is available at

table S1A. Data from cosmogenic 36Cl dating samples.

table S1B. Major element compositions of cosmogenic 36Cl dating samples.

table S1C. Trace element compositions of cosmogenic 36Cl dating samples.

table S2. Static GPS measurements of highstand paleoshoreline.

fig. S1. Histogram of elevation differences between GPS and DEM.

fig. S2. One-meter-resolution DEM and identification of highstand paleoshoreline outcrops.

fig. S3. Estimated vertical displacement in response to crustal unloading from lake draining.

fig. S4. Map of highstand paleoshoreline surface where it intersects the trace of the Troncoso fault.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank our many LdM collaborators for enlightening discussions over many years, in particular J. Fierstein and W. Hildreth. Comments by an anonymous reviewer, M. Pritchard, R. Denlinger, and C.-T. Lee helped improve this paper and are much appreciated. L. Torres, J. Torres, N. Lord, P. Sobol, L. Lara, A. Amigo, A. Alarcon, OVDAS, Servicio Nacional de Geología y Minería (SERNAGEOMIN), the U.S. Geological Survey (USGS) Volcano Disaster Assistance Program, and S. Zimmerman and A. Hidy at LLNL are thanked for logistical, analytical, and computational support. Funding: This research was supported by U.S. NSF grants EAR-1411779 and EAR-1322595, OVDAS/SERNAGEOMIN, and the USGS Volcano Disaster Assistance Program via the U.S. Agency for International Development Office of Foreign Disaster Assistance. Author contributions: B.S.S., H.L.M., and B.T. designed this study. B.S.S., H.L.M., B.T., L.C., N.L.A., N.G., and A.K.D. conducted field work including GPS measurements and photogrammetry. H.L.M., J.M.L., L.C., N.L.A., N.G., A.K.D., and K.L.F. conducted the analyses and modeling. B.S.S. and H.L.M. wrote the paper with input from all other coauthors. All authors contributed to interpretation and analysis of the data. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from B.S.S.

Stay Connected to Science Advances

Navigate This Article