Breathing of the Nevado del Ruiz volcano reservoir, Colombia, inferred from repeated seismic tomography

Nevado del Ruiz volcano (NRV), Columbia, is one of the most dangerous volcanoes in the world and caused the death of 25,000 people in 1985. Using a new algorithm for repeated tomography, we have found a prominent seismic anomaly with high values of the Vp/Vs ratio at depths of 2–5 km below the surface, which is associated with a shallow magma reservoir. The amplitude and shape of this anomaly changed during the current phase of unrest which began in 2010. We interpret these changes as due to the ascent of gas bubbles through magma and to degassing of the reservoir. In 2011–2014, most of this gas escaped through permeable roof rocks, feeding surface fumarole activity and leading to a gradual decrease of the Vp/Vs ratio in the reservoir. This trend was reversed in 2015–2016 due to replenishment of the reservoir by a new batch of volatile-rich magma likely to sustain further volcanic activity. It is argued that the recurring “breathing” of the shallow reservoir is the main cause of current eruptions at NRV.

dome material with the superficial structure. Drumbeat-type seismicity is associated with this 23 phenomenon 7 . Among its most striking features is the similarity with low energy VT events, 24 occurring with comparable waveforms and energy at relatively regular time intervals. Episodes 25 of drumbeat seismicity of short duration are still detected today. 26 The observations are consistent with a simple model of sulfur degassing 23 based on (1) 27 gas accumulation beneath an impermeable cap during times of repose, and (2) periodic gas 28 release when the cap ruptures. The continued release of SO2 after each eruption suggests that the 29 cap is gradually released over several years. 30 In juvenile solid material emitted during the 1985 and 1989 eruptions of the NRV, Stix et 31 al. 23 found evidence for pre-eruptive magma emplacement at shallow levels in (1) anhydrous 32 mineral assemblages of plagioclase and pyroxene, (2) high silica contents of glasses, and (3) low 33 water contents in melt inclusions (averaging between 1.6 and 3.3 wt.%). They proposed a 34 multistage model of magma transport and degassing that involves alternating periods of magma 35 ascent and magma ponding. Initially, volatile-bearing magma ascends from a deep reservoir 36 located at depths of 9 -15 km, driven by buoyancy. During decompression, the magma loses 37 gas, particularly CO2 and sulfur. The magma eventually ponds at its neutral buoyancy level. 38 Observations suggest a period of magma storage at shallow depths, where gas-saturated magma 39 cools and crystallizes, thereby releasing gas. As a result, CO2 is depleted from the residual melt 40 whereas H2O and SiO2 are enriched due to fractional crystallization. H2O enrichment is also due 41 in part to increased solubility in the melt as CO2 is degassed. Lastly, the density of the magma In this study, we use the catalog of arrival times of P and S seismic waves from local 47 events recorded by permanent seismic stations at Nevado del Ruiz Volcano (NRV). There are 48 57,646 events in the initial dataset, corresponding to 429,154 P-phases and 400,969 S-phases 49 recorded by 124 seismic stations operated in the NRV area at different times. In this study, we 50 used the data from 1998 onwards. Data from earlier times are available but were not used 51 because they came from a relatively small number of stations.

52
The distributions of events for several time intervals are shown in Extended Data Figure   53 2. During the observation time, seismicity was largely concentrated beneath the summit area and 54 started to migrate laterally towards the northwest and southwest in 2012. 55 We have selected data from four different time-intervals: 1998-2010, 2011-2012, 2013-56 2014, and 2015-2016, as indicated in Figure 1B by different colors. Here, we have considered 57 the differences in the tomography results in periods 1, 3 and 4 with respect to results for the same 58 period 2. This allowed us to avoid considerable differences in data distributions.

59
The data was selected to achieve maximum similarity of ray path distributions between 60 the pairs of data subsets. The selection procedure was comprised of the following steps. (1) For 61 the two-time episodes, we select one containing fewer data than another. (2) In this dataset, we 62 range the events according to the number of picks per event from the maximum value, to the 63 minimum one. (3) For the first event in the 1 st dataset, we select all events in the 2 nd dataset 64 located at a distance less than a predefined value (0.5 km in our case). (4) For the selected events 65 in the 2 nd dataset, we count the number of common phases with the current event in the 1 st subset 66 (same stations and same types of wave, P or S). (5) We select one or several events having the 67 maximum number of common phases. (6) In the case of several events with the same maximum 68 number of common phases, we select one located at a minimum distance from the current event 69 in the 1 st dataset. (7) The events are taken into consideration if the number of common phases is 70 less than or equal to 8. The selected events are removed from both initial subsets. (8) The same 71 procedure is repeated for the next event in the 1 st subset. The numbers of the selected events and 72 time pick in all cases are shown in Extended Data Table 1. In all cases, the total number of 73 involved stations was 19.

74
The inversion was performed using a modified version of the LOTOS code 26 . The 75 procedure starts with preliminary source locations using a 1D reference model and the grid 76 search method. To speed up the preliminary location, we used straight ray paths to calculate 77 travel times. The 1D reference model was the same as in the case of locating the entire dataset 78 prior to data selection.

79
The iterative inversion procedure contains the recurrent steps of source locations in the 80 3D velocity model, matrix calculation, and the inversion. The source location procedure, in this 81 case, uses the 3D bending algorithm for ray tracing 27 .

82
The 3D distributions of the P-and S-wave velocities are parametrized by a set of nodes 83 distributed in the study area according to the density of rays. Between the nodes, the velocity is 84 approximated continuously using the tri-linear interpolation. The grids are constructed in the first 85 iteration; then the velocity anomalies are updated for the same nodes. To reduce the effect of the 86 grid geometry on the result, we performed the inversions for several grids with different basic 87 orientations (0°, 22°, 45°, and 66° in our case), then computed the final model as the average of 88 the resulting distributions. It is important that the grids be created for the first data subset; for the 89 second subset, they are just copied from the first case.

90
The inversion was performed using the LSQR method 28,29 . We performed simultaneous  An important part of any tomography study is synthetic modeling, which allows an 105 assessment of the spatial resolution of the results and the determination of optimal values for the 106 inversion parameters. In this study, synthetic modeling is necessary to distinguish between the 107 impact of different ray path configurations and real velocity variations at depth.

108
The synthetic data were computed by three-dimensional ray tracing through a predefined 109 synthetic model. In the LOTOS code, there are several possibilities for defining synthetic 110 anomalies, both in map view and in vertical profiles. After calculation of synthetic travel times 111 between the actual distribution of sources and receivers, the data were perturbed with random 112 noise (0.02 and 0.05 s in our case) that provides the same variance reduction as in the case of the 113 experimental data analysis. Any information about the source coordinates and origin times was 114 "forgotten". The recovery procedure began with locating the events in the initial 1D model that usually bias the sources significantly. The iterative inversion procedure and the inversion 116 parameters used for recovering the synthetic model were identical to those for the experimental 117 data.

118
To save space, we present the results of synthetic modeling for Series 1 only. For the 119 other series, results appear to be similar and even better, because of a larger amount of data.
120 Extended Data Figure S3 shows the results of the checkerboard test, in which the starting

125
for both cases show that the main pattern is recovered (Extended Data Figure S3). In the bottom 126 row, we present the difference between the recovered models. Deviations in the central part of 127 the study area are less than 2% (note that the color scale interval for the difference is half that for 128 the recovered models). A local instability in the S-velocity model is observed in the southwestern 129 edge with a magnitude of ~4%, which we attribute to a few event mislocations and differences in 130 the ray configurations in this part of the study area. Nevertheless, the changes in the model are 131 significantly smaller than the amplitude of the main anomalies.

132
Another series of tests is presented in Extended Data Figure S4. In this case, we created 133 synthetic models with realistic anomaly distributions in a vertical section, i.e. which are similar 134 to those obtained using the real data. The anomalies were defined as polygons with unchanged 135 shapes in the direction across the section. We defined the amplitudes of the P and S anomalies Based on the experimental data, we computed a total of six models, two for each series.

152
In all cases, we performed three iterations and used the same inversion parameters. The

173
In the vertical section, we see that beneath the volcano summit, the higher P velocity 174 coexists with a strong negative S anomaly that results in a very high Vp/Vs ratio. In the first time 175 interval, it exceeds 2.2, and in the second interval, it reaches 2.0. It should be noted that this 176 anomaly of high Vp/Vs ratio matches a narrow, nearly vertical zone of seismicity just beneath the 177 summit. The difference between the models in the first series shows that the structure is mostly 178 changed by an increase of the S velocity to more than 10% beneath the summit at 2-3 km asl.

179
This causes the corresponding decrease in the Vp/Vs ratio of more than 0.3.

180
It is important to compare the inversion results for the same interval of 2011-2012 181 derived for the three series. We see some differences that are especially prominent for the S-182 velocity distribution. These differences are merely due to changes in the data configurations.

183
This indicates that comparing results without constructing identical datasets might be risky.
Interestingly, despite considerable differences in the distributions of the P and S anomalies, the 185 models of Vp/Vs ratios look more similar in this case.

186
For the second series (Extended Data Figures S7 and S8), we observe a further decrease 187 in the Vp/Vs ratio beneath the volcano. In the period 2013-2014, at 2-3 km asl, where we expect 188 the magma reservoir to be located, the Vp/Vs ratio appears to be close to the average value in the  191 In the vertical sections corresponding to all time intervals, we observe a strong shallow 192 anomaly of high Vp/Vs ratio in the summit area. This can be interpreted as strongly fractured, 193 highly saturated rocks, which remain almost unchanged during the entire observation period.

194
In the case of both series, we observe considerable changes in the P and S anomalies and  Numerical estimates for the magma sources 207 NRV is known for its very large venting of SO2 gas and provides an excellent example of 208 the general problem of the origin of volcanic gasses. In many volcanoes, it has proven to be 209 difficult to reconcile the sulfur and magmatic budgets, such that the amount of sulfur that is lost 210 is larger than what can be accounted for by magma 20 . Magma sulfur concentrations are 211 determined from melt inclusion and glass data but, by definition, only provide values for magma 212 that has already undergone crystallization. Also, it is known that sulfur solubility is low 213 compared to that of H2O, such that it may form gas at large pressures 22 . This has led to the 214 conclusion that many magmas carry exsolved gas as they rise towards shallow storage zones 21 .

215
Magma storage prior to the eruption, therefore, allows gas to escape and to get vented at the 216 surface or through a near-surface hydrothermal system.  We can also estimate the minimum magma volume using the amount of sulfur or water 229 and assuming that it got degassed passively before the 2015 dome-building eruption. The total mass of SO2 gas erupted in the 2012-2015 period, over four years, is 7x10 6 tons. The total 231 amount of sulfur in subduction zone magmas has been reconstructed using phase equilibria 232 relationships and data on both gas and melt 20 , and is between 0.05 and 0.5 wt% with very rare 233 exceptions. Using an average value of 0.1 wt%, we can estimate the amount of melt that is 234 required to supply the mass of SO2 gas that got vented. The amount of sulfur is 3.5x10 6 tons, km. An alternative calculation can be made using SO2/H2O ratios, which are typically between 238 2% and 5% by weight 20 . At the neighboring Galeras volcano, this ratio is 3% 20 . We can, 239 therefore, convert the mass of SO2 vented into one of 2.3x10 11 kg for H20. Assuming that magma 240 entered the shallow reservoir with 3.3 wt% H20 and that it progressively degassed to a 241 concentration of 1.6 wt% 23 , we obtain a magma volume of 5x10 9 m 3 and a reservoir diameter of 242 about 2 km. These estimates depend on various ratios that are not known precisely and are only 243 accurate to within a factor of about 2. The orders of magnitude, however, are well constrained 244 and lead to reservoir volumes that are consistent with our independent seismological estimate. 245 We can relate these observations to the melt budget of the reservoir. The total volume of 246 the reservoir may change as a function of the reservoir average pressure PR through deformation 247 of the wall rocks. Changes of pressure induce changes of dissolved volatile content and density, 248 such that one can define an effective compressibility , which is much larger than that of the wall 249 rocks 25 . The pressure change can be calculated as follows 24,25 : where V is the reservoir volume, PR pressure, Qin the mass flux of magma into the 252 reservoir at density , Qout the mass flux of gas leaving the reservoir at density g. The last term 253 on the right-hand side is: where T is temperature, which is positive and represents a contribution due to 256 crystallization, which acts to increase the volatile content of the residual melt and gas mixture.

257
Assuming that there is no recharge (Qin=0), for example, this equation has a steady-state solution 258 with no change of pressure such that the amount of volatile that gets exsolved due to 259 crystallization balances the amount of gas that gets vented out of the reservoir. A solution to this 260 equation requires two closure relationships specifying how the inputs and outputs (Qin and Qout) 261 vary as a function of reservoir pressure.

262
In an open conduit system connecting the shallow reservoir to a deeper magma source 263 located at vertical distance h beneath the shallow reservoir (Extended Data Figure S11), the 264 pressure difference driving the flow of magma into the reservoir can be written as: where PS is the source pressure. This shows that magma replenishment is enhanced by 267 decreasing pressure in the shallow reservoir. It also shows that replenishment ceases when 268 pressure in the shallow reservoir reaches a value equal to (PS - g h).

269
The driving pressure for gas venting, assuming permeable roof rocks, will always be 270 positive and large owing to the small density of gas. Thus, the control on the rate of degassing is 271 mostly due to the permeability, which is unknown, and to the availability of gas in the reservoir.  Degassing of the reservoir is fed by gas bubbles rising through the reservoir. With time, 287 as stated above, one expects that a degassed region grows at the base of the reservoir. If one can 288 track the rise of the boundary between the gas-depleted lower region and the gas-rich upper 289 region, one can estimate an average ascent rate for gas bubbles. According to our tomographic 290 images, this rate is about 1 km per 4 years or about 10 -5 ms -1 . Using the well-known formula for 291 the velocity of gas bubbles through melt: where a is the bubble radius,  the melt density,  the melt viscosity and where we have 294 assumed that the density and viscosity of gas are very much smaller than those of melt. For Extended Data Figure S2. The distributions of events in different times of observations. The time of the interval and the number of events are indicated above each plot. The red dots are the earthquakes, and the blue triangles are the seismic stations. Contour lines depict the relief 30 . This picture is produced using Surfer 12, Golden Software 31 . Figure S3. An example of the checkerboard test for the Series 1. Upper two rows show the inversion results at the altitude of 3 km for the P-and S anomalies and Vp/Vs ratio. The lower row is the difference between these models. The shapes of the initial synthetic anomalies are depicted with dotted lines. The contour lines depict the relief 30 . Note that the color scale intervals for the anomalies and differences are different. This picture is produced using Surfer 12, Golden Software 31 . Figure S4. Results of modeling the repeated tomography for the Series 1 with synthetic anomalies of realistic shape. Left column presents the shapes of the initial anomalies of Vp/Vs ratio for three different models. Number pairs indicate the values of P and S anomalies in percent. Middle row is the recovering results in the case of the same synthetic Model 1 and difference (lower plot). Right column presents the reconstruction results and the difference between two different models 2 and 3. This picture is produced using Surfer 12, Golden Software 31 . Figure S5. Results of experimental data inversions for two-time intervals corresponding to the series 1 at the altitude of 3 km asl. The distributions of the P and S anomalies and Vp/Vs ratios are presented. The line indicates the location of the vertical section used for presenting the main results. The contour lines depict the relief 30 . Note that the color scale intervals for the anomalies and differences are different. This picture is produced using Surfer 12, Golden Software 31 . Figure S6. Results of experimental data inversions for two-time intervals corresponding to the series 1 in the vertical section. The distributions of the P and S anomalies and Vp/Vs ratios and their differences are presented. The dots depict the earthquakes located at distances of less than 0.4 km from the profile. This picture is produced using Surfer 12, Golden Software 31 .