Analysis of airborne and terrestrial lidar data demonstrates that >0.4 km3 of magma cooled in sills at shallow (< 1 km) depth in the now-eroded Pliocene San Rafael Swell distributed volcanic field, Utah (USA). The volumes of each of seven sills are estimated from three-dimensional (3-D) models of the lidar data and range from 10−4 to 10−1 km3. Directions of magma flow during emplacement are interpreted from precise sill thickness measurements and measurements of linear vertical offsets within the sills, helping to identify feeder conduits and dikes; 3-D map relationships derived from lidar data demonstrate that magma flowed into and out of sills from these active dikes and eruptive conduits. Mapped sill volumes account for >92% of intrusive material within the 50 km2 study area. We conclude that sills played a significant role in modifying eruption dynamics during activity in San Rafael, and suggest that monitoring of sill inflation and deflation in active distributed volcanic fields may provide key information about unrest and potential eruption dynamics.


Intermediate to shallow crustal storage of pre- and syn-eruption magma modulates magma supply rate in many volcanic systems. At Mount St. Helens (A.D. 1980 eruption; USA) and Parícutin (A.D. 1943 eruption; Mexico), magma supply rate is thought to have been influenced by the presence of shallow (<10 km) temporary magma storage (Cashman and McConnell, 2005) and by the length of storage time (Scandone et al., 2007). Erlund et al. (2010) identified increasing amounts of shallow crust (≤4 km depth) assimilated at Parícutin over the 9 yr eruption, and concluded that a shallow intrusion network formed early and caused later eruptive products to be more effusive. On shorter time scales, gas may segregate preferentially into conduits above shallow sills, increasing volumetric flow in the conduit and intensifying the eruption (Conte, 2000; Pioli et al., 2009). Sill-like intrusions into shallow magma chambers have recently been geodetically linked with interferometric synthetic aperture radar (InSAR) and seismic monitoring to eruptions at Tungurahua, Ecuador (Biggs et al., 2010), and Eyjafjallajökull, Iceland (Tarasewicz et al., 2012). These models and observations suggest that it is critical to understand the volume, depth, and distribution of sills in order to forecast eruption dynamics and the evolution of volcanic systems. In young volcanic fields, such as the one around Parícutin (Connor, 1990), it is not possible to directly observe the shallow plumbing system. Here, we use lidar technology to map part of the eroded San Rafael Swell (Utah, USA) volcanic field. These data demonstrate that sills are prevalent at shallow depths (<1 km), modulated magma flow in eruptive conduits, and likely influenced eruption dynamics within this volcanic field.


The San Rafael volcanic field was active between 4.6 and 3.8 Ma (Delaney and Gartner, 1997). This volcanic field is part of a larger occurrence of Cenozoic basaltic volcanism in the Colorado Plateau and Basin and Range provinces but is distinct from many other fields as it has been eroded to a depth of ∼800 m based on its age and late Cenozoic erosion rates (e.g., Pederson et al., 2002). The sill-and-dike swarm, or volcanic plumbing system, cut a Jurassic sedimentary section from the Carmel Formation through the Cutler Formation. Diabasic dikes in this area trend 335° to 0° (relative to north) along regional joint sets, indicating low horizontal deviatoric stress during emplacement (Delaney and Gartner, 1997). Sills in the San Rafael Swell range from <5 m to >40 m thick and are exposed in cliff sides and canyons in outcrops that extend for hundreds to thousands of meters. This shallow magma plumbing system has been mapped (Delaney and Gartner, 1997) and used to improve our understanding of (1) dike emplacement (Delaney et al., 1986), (2) magma diapirism in the shallow crust (Diez et al., 2009), and (3) the spatial relationships between dikes and conduits (Kiyosugi et al., 2012), the latter of which are commonly surrounded by brecciated country rock, indicating conduit erosion during rapid magma ascent. Although Gartner (1986) described physical characteristics of exposed sills in the area, the complex emplacement processes of sills in this volcanic field have remained enigmatic. With the aid of lidar, we are able to document the complex map relationships between intrusions in the area.


Terrestrial laser scanning (TLS) performed in 2010 and 2012 collected a 7.3 GB point cloud over an area of 5 km2. Both TLS surveys used Riegl terrestrial scanners. An airborne laser scanning (ALS) survey in 2013 provided data for a 54 km2 airborne laser swath map (ALSM) (Richardson, 2013), connecting the two TLS surveys into a single study area (Fig. 1). Instrument specifications and data formats from these surveys are outlined in Table DR1 in the GSA Data Repository1. Co-registering the coordinate systems of the three surveys creates a coverage area of ∼50 km2, within which relief varies by up to 500 m. Thus, we are able to characterize the magma plumbing system in a tabular block approaching 25 km3 in volume, which bounds the study area extent. This reconstruction of a magmatic plumbing system attempts to model the amount of magma emplaced into the crust due to Pliocene volcanism for the 25 km3 space.

The three point clouds (two TLS and one ALSM) were consolidated and analyzed using LiDAR Viewer software (Kreylos et al., 2008). Because the three-dimensional point cloud is so precise, lidar data help identify subtle changes in sill thickness over large areas, vertical offsets in sills, and disrupted stratigraphy in overlying sedimentary units, which allow magma movements to be deduced. Contacts between igneous and sedimentary rocks are identified by shade contrast (igneous rocks are generally darker than sedimentary rocks in the near-infrared point cloud) and weathering patterns easily observable in the point cloud (Fig. DR2 in the Data Repository). Thickness measurements are made in LiDAR Viewer where sill upper and lower contacts are seen in close proximity. The exact locations of sill contacts are manually picked between points in the point cloud, where one point is interpreted as sill and the other as sedimentary rock (Table DR2). Uncertainty at each measurement is determined as the average of point-to-point distances on top and bottom of the sill and is drastically reduced in areas where both TLS and ALSM data are available. Other measurements made in LiDAR Viewer include sill base elevations and strikes and dips of continuous sill segments and of sedimentary host rock below sills. Locations where sills abruptly change stratigraphic level are also mapped in the field. These abrupt changes can be traced between outcrops with point cloud measurements.

Sill exposures are mapped using 1 m National Agriculture Imagery Program (NAIP) images and the ALSM digital elevation model (DEM). These are combined with thickness measurements to estimate terminal boundaries of sills. The lateral edges of sills are not commonly preserved in outcrop, so we have estimated the terminal boundary of each sill to extend no more than 0.5 km from current exposures. Sills commonly crop out at cliffs with little horizontal exposure area (Fig. 2), and by assuming that the mapped sills are contiguous between these outcrops, mapped sill areas are relatively small in comparison to interpreted sill areas (Table 1). Sill volume and average thickness are modeled by constraining the thicknesses of sills at their respective modeled boundaries to be 0 m and interpolating a Laplacian-spline surface within sill boundaries, calibrated to the measured thicknesses (Fig. 3). Results from mapped and modeled areas, and maximum measured thicknesses, are detailed in Table 1.

Using aerial images and field mapping, Kiyosugi et al. (2012) mapped 16 conduits and ∼180 vertical, en echelon dike segments, with a cumulative length of ∼53 km, that crop out in the study area (Fig. 1). The cumulative volume of igneous material stored in dikes is estimated to be the product of dike length, the modeled block height, and 85 cm, the modal dike thickness (Delaney and Gartner, 1997). This might be a slight overestimate as some dikes might not have cut through the entire block height. The volume stored in conduits is the product of the surface area, mapped with the ALSM DEM and NAIP images, of each conduit and the modeled block height. This assumes conduit thickness does not change within the vertical limits of this reconstruction and might underestimate volume if conduits widen toward the surface or formed above the present-day surface.


Seven isolated sills crop out within the study area and are vertically separated from each other by dozens of meters of sedimentary rock or horizontally by several kilometers. We interpret these sills to have been emplaced independently as a result of single dike injections, based on evidence described below. Sill volumes range from 10−4 to 10−1 km3 and have been emplaced over areas of 10−1 km2 to tens of square kilometers (Table 1). Through modeling sill geometries, we find that ∼0.4 km3 of igneous material is permanently stored in the sills, representing 93% of all intrusive rocks in our reconstructed volume. Table 2 summarizes mapped areas and modeled volumes of sills, dikes, and conduits within the study area. By combining adjacent conduits along the same dike, we estimate that 12 distinct volcanic events are represented within the study area. Emplacement processes of sills and their role in the development of the Pliocene volcanic field can be further understood by investigating individual sills.

Hebes Sill

The sill at Hebes Mountain (Fig. 1) is primarily preserved as a single 1.9 km2 sill exposed over an area of ∼4 km2. This sill generally dips with strata 1°–8° from horizontal, to the northwest, although locally some areas dip 5°–30° toward the center of the sill. Sill thicknesses are measured to a precision of ±75 cm, with virtually all exposures measuring >19 m. The sill thins monotonically from the center to the edges of Hebes Mountain, thinning most rapidly to the southwest.

By modeling Hebes sill as a 5.4 km2 area, roughly following the shape of Hebes Mountain, the volume is estimated to be 8.5 × 10−2 km3. The elongate nature of this sill (Fig. 3), with increased thickness trending in the northwest dip direction, is aligned in the regional dike direction, perhaps indicating a linear source region (dike) feeding the sill.

Central Cedar Sill

The Central Cedar sill caps two buttes to the east of Cedar Mountain and is exposed on the east-facing cliffs of Cedar Mountain (Fig. 2; Fig. DR2). The outcrops are interpreted to be parts of a single sill, as their basal contacts project across the small valleys between each exposure at the same elevation. Measured thicknesses of the Central Cedar sill are 2–15 m, with basal contacts that dip with sedimentary host rock at 2°–5° WNW to SW. The sill crops out adjacent to a conduit associated with a ∼2-km-long dike on Central Cedar Mountain. Basalt between the conduit and sill appears continuous, with no brecciation, suggesting that the dike and sill were formed cotemporally and were thus comagmatic.

The average uncertainty in thickness measurements for Central Cedar sill is <20 cm, due to coverage from both ALSM and TLS data sets. The point cloud also enables the mapping of curvilinear “step-up” features, defined by Gartner (1986) as vertical offsets between different intrusion pathways. Flow direction during intrusion is interpreted to be parallel to step-ups, indicating flow to the W–WNW, away from and/or toward, the conduit. Modeled as a tongue-shaped body intruding to the west from the suspected source dike (Fig. 3), Central Cedar sill has an areal extent of ∼0.88 km2 and a total volume of 4.4 × 10−3 km3 (Table 1).

A linearly thinning trend away from the conduit is evident in the sill, continuing for 1 km to the observed sill limit (Fig. 4). Within 100 m of the conduit, sill thickness changes dramatically due to the presence of rotated sandstone blocks with thin basalt lenses injected over the tops of the sandstone blocks, indicating roof collapse into the sill (Fig. 2). From these observations we conclude that the Central Cedar sill was fed from a single dike and was emplaced in a tongue-like fashion to the west in its initial dipping direction. Further, we infer that a conduit-forming volcanic event may have halted further advance of the sill, and subsequent flow of magma from the sill into the conduit caused the observed conduit-adjacent roof collapse.


Through lidar mapping of the San Rafael study area, seven sill-forming events in the shallow crust and 12 conduit-forming events are identified and mapped in detail (Fig. 1). We model the total volume of igneous material stored in sills to be 0.4 km3 within a 25 km3 block. This sill volume represents 93% of the stored igneous volume in the block, with the remaining 7% in dikes and conduits. There is no doubt that, volumetrically, sills are a critical component of the magma plumbing system in this distributed volcanic field.

It is possible, in fact, that sill volume in the San Rafael Swell volcanic field is comparable to erupted volume. Eruption volumes for the 12 conduits cannot be directly observed, as those lavas are completely eroded away. Eruption volumes for monogenetic volcanoes in similar fields span three orders of magnitude, ranging from 10−3 to 10−1 km3 (e.g., Crowe et al., 1983; Condit et al., 1989). If we assume that average eruption volume is 0.1 km3 for conduits in the San Rafael Swell, ∼1.2 km3 of basalt would have been erupted at the surface, four times the estimated sill volume. Again, this comparison suggests that, volumetrically, crustal storage of magma in sills is a major feature of the magmatic system.

The general shape of the mapped sills in this field does include a thick center of several to tens of meters in height, tapering edges, and horizontal dimensions of one to several kilometers. While smaller sills exhibit a monotonic decrease in thickness from their interiors (Fig. 4), the thickness profiles of larger sills are more complex and multiple thickened zones exist (Fig. DR1). The irregular shapes and size range of these sills might suggest that all sills in this area are the product of single injection events and are not polygenetic (Gudmundsson, 2012). Furthermore, the maximum observed thicknesses of each of these sills are not highly correlated to the exposed or modeled areas of each sill. Sills in this area generally ascend stratigraphy only after lifting the roof, enabling exploitation of new bedding planes, suggesting that initial emplacement at this depth was a pressure-driven process, similar to the intrusion of the Trachyte Mesa laccolith (Henry Mountains, Utah, USA) (Wetmore et al., 2009). Because horizontal deviatoric stress was low in this area during the Pliocene (Delaney et al., 1986), the minimum compressive stress direction could have significantly migrated from horizontal at 1 km depth, enabling sill formation given local stress conditions influenced by overlying topography (Gudmundsson, 2012).

The development of shallow sills likely affected eruption dynamics. If a comagmatic sill is present during a volcanic eruption, ascending bubbles can become concentrated in the vertical conduit at the conduit-sill junction by disproportional capture of the liquid phase of a two-phase flow in the horizontal branch (Conte, 2000). This concentration occurs if overall magma flux is sufficiently low. The presence of the sill, therefore, enables modulation of explosive potential. Following the method of Pioli et al. (2009), assuming a magma density of 2800 kg/m3, we calculate the transition flux to be 1.8 × 105 kg/s within the Central Cedar Mountain conduit (diameter 25 m), where lower flux would have concentrated bubbles in the conduit system. As average mass eruption rate for Strombolian eruptions is commonly observed to be 103–105 kg/s (Pioli et al., 2009), the presence of sills at this level, where H2O exolves, critically impacts volcanic hazard.

Sills have been identified as a major instigator of unrest in association with stratovolcanoes (e.g., Biggs et al., 2010; Tarasewicz et al., 2012), volcanic calderas (Macedonio et al., 2014), and monogenetic volcanic eruptions (Erlund et al., 2010). The observation in the San Rafael Swell that the vast majority of igneous rock at ∼1 km depth is contained in sills suggests that similar deformation events may be precursory to volcanic eruptions in some active volcanic fields. Monitoring of active volcanic fields may benefit from use of deformation networks to detect sill emplacement in the shallow crust.

This project would not have been possible without the support of D. Phillips and UNAVCO for lidar field equipment and data processing assistance. TLS data acquisition was funded by a grant from the National Science Foundation (NSF) (EAR-0910696). ALS data acquisition and processing were completed by the National Center for Airborne Laser Mapping (NCALM) with funding provided by NSF’s Division of Earth Sciences, Instrumentation and Facilities Program (EAR-1043051). M. Oskin and O. Kreylos of University of California–Davis are thanked for their help regarding LiDAR Viewer. University of South Florida students and faculty J. McIlrath, K. Mannen, K. Kiyosugi, J. Wilson, S. Tavarez, and T. Doering are thanked for field assistance. Helpful comments from K. Cashman, A. Gudmundsson, and an anonymous reviewer improved this manuscript.

1GSA Data Repository item 2015345, Table DR1 (lidar equipment specifications), Table DR2 (sill thickness measurements made in LiDAR Viewer), Figure DR1 (maps of additional sill thickness models), and Figure DR2 (photograph and point cloud perspective of igneous bodies in study area), is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.