Causes of volcanic unrest at Mt. Spurr in 2004–2005 inferred from repeated tomography

Mt. Spurr is the largest active volcano in Alaska of high explosive potential. The most recent activity, including two recent magmatic eruptions in 1953 and 1992, has occurred via the flanking Crater Peak. From 2004 to 2006, strong seismicity, gas flux, and heating were observed in the summit area, which had remained inactive for more than 5 Ka. To understand the cause of this reactivation, we performed repeated tomography inversions that clearly imaged the magma reservoir beneath Mt. Spurr and showed temporal changes in its shape and intensity. During the two years preceding the unrest, we observed ascension of the upper limit of the reservoir-related anomaly from a depth of 5 to 3 km below the surface, accompanied by strong seismicity. During the following years, the shape of the anomaly remained unchanged, but its intensity weakened. These observations may indicate the release of fluids from the ductile reservoir and fast upward ascent through the brittle cover that caused intensive seismicity and gas flux during the unrest from 2004 to 2006. The origin of this zone will possibly cause a resumption of explosive eruptions in the summit area of Mt. Spurr.

Details of the algorithms of repeated tomography including the data selection are 22 described in Vargas et al. 23 . For this study, we selected three time intervals: [1996][1997][1998][1999][2000][2001] (1), 23 Figure 2 of the main paper. This selection was 24 done regarding the different stages of volcano activity, but also considering the distribution of 25 events such that the following procedures provided sufficient data for tomography. We 26 performed repeated tomography inversions for two pairs of datasets, namely (1) -(2) and (2) -

27
(3). We selected pairs of events that occurred in each of the two datasets and at distances less 28 than 0.5 km. For these, we considered the picks recorded by the same stations and having the 29 same type of waves, P or S. Among all the combinations of paired events, we selected those 30 having the maximum number of common phases. In the case of several pairs with the same 31 number, we selected events at minimum distances. For tomography, we used pairs of events with 32 a number of common phases less than or equal to 6. For the pair (1) -(2), we selected 524 paired 33 events with a number of common phases of 2,690 and 1,436 for the P and S waves, respectively. 34 In this case, 11 stations were involved. For the pair (2) -(3), the number of events was 1,325, and 35 the number of the P and S picks were 7,076 and 4,900, respectively. The distribution of events 36 together with the P and S ray paths for these cases are shown in Figure S1.

37
The tomography procedure is based on the modified version of the LOTOS code 29 . All 38 the controlling parameters, such as grid spacing, smoothing and amplitude damping, were 39 identical for all models. The most important parameters are presented in Table S2. 40 The workflow starts with the location using the one-dimensional (1D) velocity model.

41
During the preliminary stage, the source coordinates are determined using the grid search 42 method. To accelerate the calculations, the travel times during this stage were computed based on a straight line approximation. Note that prior to data selection, all events were relocated using 44 this algorithm.

45
The recurrent tomographic procedure starts with another stage of source locations using 46 the bending algorithm for ray tracing 30 that uses the starting 1D model during the 1 st iteration and updated three-dimensional (3D) velocity models for the following iterations.

48
The 3D distributions of the P and S velocity anomalies were parameterized using a set of 49 nodes distributed in the study volume according to the ray path coverage. The parameterization 50 grid was constructed during the first iteration and remained unchanged during the following 51 iterations. The parameterization nodes were installed according to the distribution of rays. No 52 nodes were installed in the ray density was less than 0.1 of the average value. In map view, the 53 nodes were set regularly with the spacing of 2 km. In the vertical direction, the distance between 54 nodes was inversely proportional to the ray density, but was not smaller than 1 km (see Table   55 S2). To make the solution grid independent, we performed inversions in several grids having 56 different basic orientations and then averaged the derived results. A special feature of the 57 repeated tomography is that the grid files are created for one dataset and then are copied to the 58 folder of another dataset.

59
The sensitivity matrix was computed along the ray paths derived after the source location  Table S2. Note that these 65 parameters were identical for all models presented here. Note that in this study, we used much less data than that used in Koulakov et al. 16 which determined the crustal structure beneath the 67 same volcano using the data of the temporary network. Therefore, here we used more 68 conservative damping values than those in the previous study which was important in obtaining 69 stable temporary variations.

70
After computing anomalies in all grids, we combined them using a regular mesh and then 71 used the 3D model as a basic distribution for the next iteration. In total, we used five iterations. were reduced for the P-and S-wave data from 0.13 s and 0.18 s to 0.079 s and 0.087 s (41% and 74 52%), respectively. During other periods, the variance reduction was similar or greater (see 75   Table S3).

76
All the results presented in the paper can be reproduced using the data files and the 77 program codes available at http://www.ivan-art.com/science/LOTOS/repeatomo.zip. This  Based on the experimental data, we computed a total of four models, two for each series.

83
In all cases, we performed five iterations and used the same inversion parameters. The   Figure S5 shows the synthetic test with realistic anomalies representing the general 132 patterns observed in the case of the experimental data inversion. Beneath the summit area, we 133 defined a columnar anomaly with contrasted high Vp/Vs anomaly. As in the previous series of tests, we started with recovering identical models (1 and 3 rows). As in the previous case, the 135 recovered distributions of the Vp/Vs ratio appear identical, and their difference is near zero.