Back to the deformation and Stress Change Modeling home page

Yellowstone National Park (Fig. 1) is famous for its numerous hydrothermal features and other natural wonders but is better known among earth scientists as the site of the world's largest restless caldera. Many earth scientists believe the park is the present day terminus of the active Yellowstone hotspot. The hotspot track can be traced along a string of large calderas (for example, 1) to its origin ~16 million years ago in southeastern Oregon and northern Nevada. The three most recent caldera-forming events occurred during cataclysmic rhyolite ash-flow tuff eruptions during the last 2.1 million years in the Yellowstone region (1, 2). Yellowstone caldera, the youngest of the three, (Fig. 1) formed ~630,000 years ago in an eruption (many times larger than any historic volcanic eruption) that ejected ~1,000 km³ of debris -- about 1,000 times the volume of magma erupted at Mt. St. Helens in May 1980 (1). A subsequent episode of dominantly extrusive volcanism buried Yellowstone caldera under rhyolite lava flows from 150,000 to 70,000 years ago (2). Even though the last magmatic eruption occurred ~70,000 years ago, geologic and geophysical evidence suggest that a crustal magma reservoir beneath Yellowstone is maintained in a partly molten state by episodic intrusions of basaltic magma (3). Because another caldera forming eruption is almost inevitable, though not imminent, a continuous monitoring program is important. It is equally important, however, to assess the patterns of deformation of the caldera, in an effort to understand the magmatic plumbing beneath Yellowstone caldera and large calderas in general.

The main advantage of satellite interferometric synthetic aperture radar (InSAR) is its ability, under favorable conditions, to measure deformation of the entire caldera floor and its surroundings. By measuring a surface of deformation rather than movement at isolated points, we are able to locate and characterize deformation sources better than has been possible to date. High resolution synthetic aperture radar (SAR) images of Yellowstone were formed using raw European Space Agency satellite radar data. The interferograms (Fig. 2) were made by taking the phase difference of pairs of SAR images, correcting for the topography [for example, see (12)] using a digital elevation model (DEM) (Fig. 1), then correcting for orbital errors (13). We selected radar images acquired only during the summer months to minimize the decorrelation between images resulting from snow accumulation. The summer images were then combined (Table 1) to pick paired images that minimized the effects of uncorrected topography [resulting from inaccuracies in the DEM, (14)] and the possibility of base-line decorrelation [resulting from a large separation in the orbits, (15)]. 

Fig. 2 (view larger GIF of figure 2, ~480 KB) Six interferograms constructed using the two-pass method (for example, 12, 25). The range of colors from violet to red, shown in the color bar on the top of (B), corresponds to one cycle of phase from 0 to 2p (one fringe) representing ~28 mm of displacement between a point on the ground and the satellite. The time interval spanned by each interferogram is referenced to the horizontal time axis where the beginning of each year is marked and labeled. The red, green and blue bars attached to each interferogram indicate the corresponding time interval. The first three interferograms are sequential whereas the last three are cumulative over longer time spans (with some overlap). Where coherence is adequate, the interferograms show that deformation extends only slightly beyond the northeast caldera boundary to the northwest-southeast trending faults mapped by Christiansen (Fig. 1, (2)). (A) August 1992 to June 1993. This image shows over 30 mm of inferred subsidence centered in the northeast half of the caldera and closely associated with the SC dome. (B) June 1993 to August 1995. In this image the center of deformation has shifted to the southwest half of the caldera with with over 40 mm of subsidence associated with the ML dome. (C) August 1995 to September 1996. The main deformation mode is now uplift (~20 mm) associated with the SC dome [note the reversed color sequence toward the center of SC dome, relative to (A) and (B)]. The ML dome still appears to be subsiding slightly during this time interval. (D) August 1992 to June 1995: This image shows ~60 mm of subsidence. (E) July 1992 to August 1995 also showing ~60 mm of subsidence. (F) July 1995 to June 1997: This image shows over 30 mm of uplift.

To estimate the depth of the deforming body responsible for the observed subsidence, we model the interferogram in Fig. 2B with two tabular sill-like bodies (16). We chose to model Fig. 2B because it has the best combination of a long time interval and good overall image coherence. We also chose to use two sill-like bodies, because the pattern of fringes in Fig. 2B requires more than one source. The depth of the larger body in the model (Fig. 3A, better constrained than the small body) is 8.5 ± 4 km (95% confidence intervals). The depth coincides with that at which Miller and Smith (17) find low seismic velocity bodies beneath both resurgent domes. Previous studies placed the body at a depth of 10 ± 5 km (8) and 3-6 km [possibly extending to 9 km beneath ML (18)]. The model requires a decrease in volume of 0.035 to 0.057 km³ from June 1993 to August 1995, or a volumetric decrease rate of 0.016 to 0.027 km³/yr. These rates are comparable in magnitude to volumetric increase rates of 0.01 to 0.028 km³/yr, estimated to explain uplift of the caldera from 1923 to 1984 (18).

We suggest that interacting fluid reservoirs exist beneath the SC and ML domes, with inflation of each being regulated by flow through two conduits, one beneath SC (conduit 1) and the other connecting SC and ML (conduit 2). The opening of each conduit requires a critical pressure threshold (P_C), perhaps the pressure, P, required to overcome the lithostatic load that would tend to close the conduit walls. Inflation initiates beneath SC when influx of fluids from depth causes P to exceed P_C at conduit 1, with surface uplift occurring at the SC dome. Fluid flow increases pressure in the SC reservoir until it reaches P_C at conduit 2, initiating flow into the ML reservoir and surface uplift above it. Since P=P_C in the SC reservoir, uplift ceases there while the ML reservoir inflates. When P reaches P_C in both reservoirs, and if flow through conduit 1 continues, both should continue inflating. However, if fluid supply rate through conduit 1 decreased, inflation and concomitant uplift would decline accordingly.

The direction of subsidence migration across the caldera is not known (24) and an outlet during subsidence at conduit 1 would require that the density contrast between fluid and wallrock change sign (22). So, we suggest a third conduit as an outlet from the caldera, feeding the shallow hydrothermal features in the northwest corner of the park. In two interferograms (Figs. 2, B & F) the fringes cross the northwest caldera rim suggesting a connection between the caldera system and the northwest corner of the park. Moreover, both the change from uplift to subsidence in 1985 and the change from subsidence to uplift in 1995 are either initiated or accompanied by the two largest earthquake swarms ever recorded in the park (since 1972, in the USGS-Univ. of Utah earthquake catalog). The two swarms were located just outside of the northwest caldera rim.

REFERENCES AND NOTES

1. R. B. Smith and L. W. Braille, J. Volc. Geotherm. Res. 61, 121 (1994).
2. R. L. Christiansen, in Explosive Volcanism: Inception, Evolution, and Hazards (National Academy, Washington D. C., 1984), pp. 84-95.
3. R. O. Fournier and A. M Pitt, Trans. Geothermal Resources Council 1985 International Symposium On Geothermal Energy, C. Stone, Ed. (Geothermal Resources Council, Honolulu, HI, 1985), pp. 319-327
4. W. W. Locke and G. A. Meyer, J. Geophys. Res. 99, 20079 (1994).
5. K. L. Pierce, K. P. Cannon, G. A. Meyer, Eos 78, 802 (1997).
6. J. R. Pelton and R. B. Smith, Science 206, 1179 (1979); J. Geophys. Res. 87, 2745 (1982).
7. D. Dzurisin and K. M. Yamashita, J. Geophys. Res. 92, 13753 (1987).
8. D. Dzurisin, J. C. Savage, R. O. Fournier, Bull. Volc. 52, 247 (1990).
9. D. Dzurisin, K. M. Yamashita, J. W. Kleinman, ibid. 56, 261 (1994).
10. F. Arnet et al., Geophys. Res. Lett. 24, 2741 (1997).
11. C. M. Meertens and R. B. Smith, ibid. 18, 1763 (1991).
12. D. Massonnet, P. Briole, A. Arnaud, Nature 375, 567 (1995).
13. Fringes from orbital errors (inaccuracies in satellite positions) have been corrected for by subtracting a plane from each interferogram. The dominant orbital errors appear as nearly parallel fringes, and where coherence is good across the image the plane is easily found by inspection. We used inspection to set bounds for a grid search, then accepted the plane that gave a minimum variance-corrected interferogram. Because of poor coherence, the planes used to correct images in Figs. 2, E & F were more qualitative in nature (there were no minimum variance solutions). The magnitude of tilt across each image is: 30 mm (Fig. 2A), 300 mm (Fig. 2B), 220 mm (Fig. 2C), 300 mm (Fig. 2D), 220 mm (Fig. 2E) and 110 mm (Fig. 2F). These magnitudes are much larger than any signature that could be expected from regional tilting.
14. We assessed the DEM inaccuracies by forming interferograms from SAR images paired from summers of the same year, where the pairs had low altitudes of ambiguity (~25 m). These pairs were sensitive to DEM inaccuracies and (because of the short time interval) fairly coherent. Inaccuracies appear to be generally less than 10 m with the exception of a small area off the southeast tip of Hebgen Lake (HL) that is in error by ~25 m. This DEM artifact has the largest effect in Fig. 2B (Table 1) where it is appears as an area with half a fringe of signal which could be interpreted as ~14 mm of line-of-sight subsidence. The effect of the DEM artifact near HL is smallest in Fig. 2A (Table 1) where it would appear as an area with only ~1.5 mm of subsidence.
15. H. A. Zebker and J. Villasenor, IEEE Trans. Geosci. Rem. Sensing 30, 950 (1992).
16. An initial model was obtained by forward modeling to fit ~400,000 data pixels in Fig. 2B. Seven model parameters were estimated for each of the two sills: 1. East location of sill origin. 2. North location of sill origin, (the origin of the sill specifies the geographic position of the sill by specifying one corner of a rectangle). 3. Sill dip. 4. Sill strike. 5. Sill length. 6. Sill width. and 7. Amount of compaction (directed perpendicular to the sill surface). A Monte Carlo search was first used to refine the forward model, then used to estimate bounds on the model parameters. Any model with a variance increase (relative to the best-fitting model) that an F test revealed to be significant at the 95% level was taken as an acceptable model -- indistinguishable from the minimum variance solution at the 95% level. The range of model parameters spanned by the acceptable models was used to constrain the depths and volume changes given in the text. The dip of the sill-like bodies is 0 ±12°
17. D. S. Miller and R. B. Smith, in preparation.
18. D. W. Vasco, R. B. Smith, C. L. Taylor, J. Geophys. Res. 95, 19839 (1990).
19. Because the interferograms in Fig. 2 have varying amounts of incoherence (indicated by the amount of image speckle) on top of a long wavelength signal, we averaged the interferogram phases (which are all measured modulo 2p) over a number of pixels. In most of the caldera area, this enables us to unwrap the phases. Unwrapping translates the interferogram phases into line of sight displacement by taking care of the modulo 2p uncertainty in the phases. We averaged the phase values over a square area of pixels until the unwrapped phases of the pixels occupied by the benchmarks is similar to the unwrapped phases of each of the surrounding eight pixels. We then increased the number of pixels to the next odd number squared and repeated the phase unwrapping. For each line in Fig. 4 labeled InSAR, we plotted displacement from 36 displacement estimates (the benchmark pixel and the surrounding 8 pixels with 4 different pixel averages). For Fig. 4, A through C, the uncertainties inferred from this procedure (the spread of the lines) is largest from 0 to 12 km where the leveling line runs along the edge of YL (the lake produces essentially random interferogram phases). The heavily speckled nature of Fig. 2F from 0 to 12 km along the leveling line appears to be the cause of the poor estimates of displacement there.
20. In comparing leveling and InSAR displacements we are assuming the satellite line-of-sight measurements in Fig. 4 represent predominantly vertical deformation. However, the satellite line-of-sight is inclined 23° to the vertical, so horizontal displacement may be mapped into range displacement. We assessed the magnitude of the error introduced in using this assumption by computing the effect of horizontal displacements in the synthetic (model) interferogram (Fig. 3B). Results suggest that the largest error along the leveling route introduced by this assumption is less than 3 mm, located where the route crosses the edges of the smaller sill-like body (Figs. 1 and 2).
21. H. Tarayre and D. Massonnet, Geophys. Res. Lett. 23, 989 (1996);H. A. Zebker, P. A. Rosen, S. Hensley, J. Geophys. Res. 102, 7547 (1997).
22. The volume flux rate (V) of laminar Poiseuille flow in a cylindrical pipe is given by
    
    where µ is fluid viscosity, r s is the density of the solid wallrock,r _l the fluid density, g the acceleration of gravity and R the equivalent pipe radius [see, D. L Turcotte and G. Schubert, Geodynamics: Application of Continuum Physics to Geological Problems, (John Wiley and Sons, New York, 1982) p. 240]. Our volume flux rate is about 0.016 km³/year, the driving buoyant pressure gradient (r s - r _l)g=2000 Pa/m, and m is 10⁸ Pa-s (rhyolite), 1 Pa-s (basalt), and 10-4-10-6 Pa-s (water and vapor respectively).
23. We use the term magmatic fluids to refer to magma (molten rock) or volatile gases and liquids released from magma, whereas we use the term hydrothermal fluids to refer to steam and hot water mostly of meteoric origin.
24. From the snapshots of subsidence (Fig. 2) we cannot discern whether the subsidence deformation front moves from southwest to northeast or northeast to southwest. Because the inception of uplift was captured in Fig. 2C, however, we can conclude that the uplift deformation front moves from the northeast to the southwest (Fig. 2F).
25. D. Massonnet and T. Rabaute, IEEE Trans. Geosci. Rem. Sensing 31, 455 (1993).
26. K. Feigl and E. Dupré, Comp. and Geosci., in press.
27. We thank our colleagues R. L. Christiansen and J. C. Savage and two anonymous reviewers for thorough technical reviews of the manuscript. Helpful discussion with R. O. Fournier, D. J. Stevenson, E. Brodsky, D. Massonnet and N. Pourthie is acknowledged. We also thank P. Wessel and W. Smith for GMT software [EOS 72, 441 (1991)]. This research was supported by the Volcano Hazards Program, U. S. Geological Survey. One of us (W. T.) was supported by a visiting professorship at the Division of Geological Sciences, California Institute of Technology while some of the work reported here was carried out.

  Table 1. Information on the individual Satellite images used to construct the six interferograms shown in this paper. The phase differences are relative to the Master image. The altitude of ambiguity h_a (25), determined by the separation in orbits, is the amount of topographic change required to generate one interferometric fringe.