# The 2020 Eruption and Large Lateral Dike Emplacement at Taal Volcano, Philippines: Insights From Satellite Radar Data

M G Bato, P. Lundgren, V. Pinel, R. Solidum, A. Daag, M. Cahulogan

► To cite this version:

M G Bato, P. Lundgren, V. Pinel, R. Solidum, A. Daag, et al.. The 2020 Eruption and Large Lateral Dike Emplacement at Taal Volcano, Philippines: Insights From Satellite Radar Data. Geophysical Research Letters, 2021, 48 (7), 10.1029/2021GL092803. ird-04996271

HAL Id: ird-04996271

https://ird.hal.science/ird-04996271v1

Submitted on 18 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters

## RESEARCH LETTER

10.1029/2021GL092803

**Key Points:**

* We present a comprehensive interferometric synthetic aperture radar (InSAR)-based data, analyses, and models of Taal's pre- to post-eruptive state

* During the eruptive crisis, Taal's magma reservoir lost $0.531 \pm 0.004 \text{ km}^3$ of volume while a $0.643 \pm 0.001 \text{ km}^3$ lateral dike was emplaced

* Low-latency InSAR-derived products provided crucial and significant information to PHIVOLCS during the January 2020 eruptive event

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. G. Bato,
bato@jpl.nasa.gov

**Citation:**

Bato, M. G., Lundgren, P., Pinel, V., Solidum, R., Daag, A., & Cahulogan, M. (2021). The 2020 eruption and large lateral dike emplacement at Taal volcano, Philippines: Insights from satellite radar data. *Geophysical Research Letters, 48*, e2021GL092803. https://doi.org/10.1029/2021GL092803

Received 4 FEB 2021
Accepted 8 MAR 2021

© 2021. Jet Propulsion Laboratory, California Institute of Technology. Government sponsorship acknowledged.

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.

<figure><img src="image_2.png" /></figure>

## The 2020 Eruption and Large Lateral Dike Emplacement at Taal Volcano, Philippines: Insights From Satellite Radar Data

M. G. Bato¹ ID, P. Lundgren¹ ID, V. Pinel² ID, R. Solidum Jr.³, A. Daag³, and M. Cahulogan³

¹Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, ²Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, Grenoble, France, ³Department of Science and Technology (DOST), Philippine Institute of Volcanology and Seismology (PHIVOLCS), Quezon City, Philippines

**Abstract** On January 12, 2020, Taal volcano, Philippines, erupted after 43 years of repose, affecting more than 500,000 people. Using interferometric synthetic aperture radar (InSAR) data, we present the pre- to post-eruption analyses of the deformation of Taal. We find that: (1) prior to eruption, the volcano experienced long-term deflation followed by short-term inflation, reflecting the depressurization-pressurization of its ~5 km depth magma reservoir; (2) during the eruption, the magma reservoir lost a volume of $0.531 \pm 0.004 \text{ km}^3$ while a $0.643 \pm 0.001 \text{ km}^3$ lateral dike was emplaced; and (3) post-eruption analyses reveal that the magma reservoir started recovery approximately 3 weeks after the main eruptive phase. We propose a conceptual analysis explaining the eruption and address why, despite the large volume of magma emplaced, the dike remained at depth. We also report the unique and significant contribution of InSAR data during the peak of the crisis.

**Plain Language Summary** Taal volcano in the Philippines erupted on January 12, 2020. Here, we present the pre-, co-, and post-eruption data, model, and analyses using interferometric synthetic aperture radar (InSAR) data acquired by various satellite systems. We find that: (1) prior to the eruption, the volcano experiences a sequence of long-term (>1 year) deflation followed by short-term (≤1 year) inflation as a result of the depressurization-pressurization of its ~5 km depth magma reservoir; (2) during the eruption, the magma reservoir lost a volume of $0.531 \pm 0.004 \text{ km}^3$ while a $0.643 \pm 0.001 \text{ km}^3$ lateral dike was emplaced; and (3) post-eruption analyses reveal that the magma reservoir is in recovery starting ~3 weeks after the main eruptive phase. We propose a conceptual analysis to explain the 2020 Taal eruption and the dike emplacement. We also report the unique and significant contribution of remote sensing data, particularly InSAR during the peak of the crisis.

## 1. Introduction

Taal volcano is a caldera located in southwestern Luzon island in the Philippines (Figure 1). It is one of the frequently erupting volcanoes in the country, having at least 33 known historical eruptions between CE1572 and CE1977 (Delos Reyes et al., 2018). The volcano is part of the Macolod Corridor, a complex NE-SW trending 50-60 km-wide rift zone, characterized by active volcanism, crustal thinning, extensive faulting, and block rotation (Delos Reyes et al., 2018; Galgana et al., 2014; Pubellier et al., 2000; Wolfe & Self, 1983). Recent activities at Taal are limited to Volcano Island, a 5 km wide, 311 m high resurgent dome that has >40 volcanic vents, with the largest being the Main Crater (Delos Reyes et al., 2018; Zlotnicki et al., 2018). Geophysical studies reveal the presence of a large hydrothermal system at 2.5 km depth below the Main Crater, sitting above a magmatic body (Alanis et al., 2013; Fikos et al., 2012; Yamaya et al., 2013; Zlotnicki et al., 2018). Eruptions are controlled by local tectonics, and magma reservoir-water interactions from external sources (e.g., rainwater, groundwater, and seawater) that can result in different eruption styles and a wide range of explosivity (Delos Reyes et al., 2018; Pubellier et al., 2000; Zlotnicki et al., 2018).

On January 12, 2020, Taal volcano experienced a phreatic-phreatomagmatic eruption after 43 years of quiescence, spewing a column of steam-laden ash reaching up to ~15 km-high with frequent volcanic lightning (PHIVOLCSa, 2020). In Figure 1c, we summarized the chronology of events as the crisis progressed and waned, based on published PHIVOLCS bulletins and materials (Martinez-Villegas, 2020). A total of 565,005
<figure><img src="image_3.png" /></figure>
individuals were directly affected by this crisis with ~$69M worth of damage to infrastructure and agriculture reported (NDRRMC, 2000).

Here we present the pre- to post-eruption analyses of the January 2020 Taal eruption and dike emplacement using interferometric synthetic aperture radar (InSAR) datasets. We discuss a conceptual analysis on what drove the 2020 eruption as well as explain why the dike remained emplaced at depth. We also report the significant contribution of InSAR data during the on-going crisis when most of the in situ instruments were not operating.

## 2. Data and Methods

### 2.1. InSAR and Time-Series Processing

#### 2.1.1. Pre-Eruptive Data Set

We explored the surface deformation using synthetic aperture radar (SAR) data from two systems: the Japan Aerospace Exploration Agency (JAXA) ALOS-2 satellite's PALSAR-2 L-band (23.8 cm wavelength) sensor and the Copernicus Sentinel-1A/B dual satellites' C-SAR C-band (5.5 cm wavelength) sensor operated by the European Space Agency (ESA). All datasets were processed into differential interferograms of surface deformation projected into the radar line-of-sight (LOS) after correcting for Earth curvature and topographic effects.

The ALOS-2 SAR data acquired from December 2, 2018 to January 12, 2020 were processed at the Jet Propulsion Laboratory (JPL) using the InSAR Scientific Computing Environment (ISCE) software (Rosen et al., 2018). The Sentinel 1A/B SAR data covering the periods between October 19, 2014 and January 11, 2020 were processed through JPL's Advanced Rapid Imaging and Analysis (ARIA) Project (Bekaert et al., 2019) based on nearest three acquisitions (i.e., *N* = 3), leading to temporal baselines of up to 36 days over Taal. ARIA provides rapid, free, and open-access standard Sentinel-1 geocoded unwrapped (GUNW) interferograms that are automatically processed in the cloud using ISCE. We leverage these GUNWs since they are ready for ingestion into postprocessing algorithms (Buzzanga et al., 2020), enabling the fast generation of deformation time-series.

We computed InSAR time-series for both the ascending and descending tracks of the two satellite systems using the MintPy software (Yunjun et al., 2019). In addition, for comparison-only, a time-series was computed for the Sentinel-1 descending track (i.e., May 5, 2016–January 9, 2020) using the New Small Baseline Subset (NSBAS) processing chain, modified to allow TOPSAR data ingestion (Doin et al., 2011; Grandin, 2015). The rationale is to ensure that the neighborhood pairing (*N* = 3) used to generate the ARIA products has no effect on the resulting time-series, as short-temporal pairs can sometimes induce biases (i.e., phase misclosure, [De Zan et al., 2019]).

#### 2.1.2. Co-Eruptive Datasets

We utilized SAR data acquired by Sentinel 1A/B from 9–17 January 2020 covering the main eruptive event to produce the InSAR phase and pixel-offset maps using ISCE. The descending data set captured the deformation from 9–15 January 2020, whereas the ascending data set spans 11–17 January 2020. We derived the horizontal and vertical displacements following the method of Wright et al. (2004).

**Figure 1.** (a) Map showing Taal volcano and the local tectonic features in the vicinity including parts of the Macolod Corridor (modified after Delos Reyes et al., 2018; Pubellier et al., 2000; Rimando & Knuepfer, 2006). VI: Volcano Island; PRV: Pansipit River Valley; MCL: Main Crater Lake. White outlines are town boundaries. LM: Lemery; AG: Agoncillo; LR: Laurel; TL: Taal; SN: San Nicolas; TS: Talisay; CL: Calamba; TN: Tanauan. (b) The Philippine archipelago and the major tectonics controlling the region. Two opposing subduction zones create the oblique convergence in the Philippines: to the west are the Manila-Negros-Sulu and Cotabato trenches, and to the east is the Philippine Trench. The result of this oblique plate motion is the Philippine Fault Zone which is a 1,200-km long left-lateral strike-slip fault that traverses the archipelago from north to south (Aurelio, 2000). Orange solid lines are active faults and blue broken lines are offshore extension of the active faults. The red triangle is the location of Taal volcano. MT: Manila Trench, NT: Negros Trench, ST: Sulu Trench, ELT: East Luzon Trough, PFZ: Philippine Fault Zone, CT: Cotabato Trench, PT: Philippine Trench. (c) Summarized chronology of events based on issued PHIVOLCS bulletins and publicly available materials (Martinez-Villegas, 2020).
### 2.1.3. Post-Eruptive Datasets

We explored Sentinel-1 SAR data from January 15 to June 27, 2020 for both the descending and ascending orbits. The interferograms were processed using ISCE and we computed the post-eruptive time-series for both tracks using MintPy. We divided the time-series into two epochs (i.e., 15 January to 04 February and 02 February to 27 June) based on the change in deformation trend within Volcano Island and derived the horizontal and vertical cumulative displacement maps following Wright et al. (2004).

## 2.2. Geodetic Source Modeling

For modeling, we used different software, algorithms, and techniques depending on the data set and the process that we intend to address (i.e., pre-eruptive: inflating magmatic source, co-eruptive: deflating magma reservoir and dike emplacement, post-eruptive: further propagation of the dike). For a detailed description of our methods, please refer to the Supplementary Information.

2.2.1. Pre-Eruptive Modeling

We used the Caltech-JPL-developed AlTar v2.0 Bayesian inversion software. AlTar is based on the Cascading Adaptive Transitional Metropolis in Parallel method of Minson et al. (2013). To constrain the model, we utilized the Sentinel-1 and ALOS-2 velocity maps from the time-series inversion covering one year of the inflation event before the eruption. We subsampled each data set and computed the covariances by applying a model-based quadtree approach (Lohman & Simons, 2005) using a sill model, which is appropriate for an inflationary volcanic signal.

We used the compound dislocation model (CDM) (Nikkhoo et al., 2017) in AlTar to model the magmatic reservoir beneath Taal. The CDM approximates an arbitrarily oriented and shaped ellipsoidal reservoir (from cigar to pancake). Overall, the CDM has 10 parameters: three for the locations (x, y, z), three for the semi-axes lengths (a, b, c), three for the rotations (ωx, ωy, ωz), and a uniform opening, u. We also estimated the shifts in the InSAR data since InSAR measurements are relative measurements based on a reference point. In the end, we estimated 14 parameters. We extract the optimal parameter values by calculating the maximum *a posteriori* probability (MAP) solutions.

2.2.2. Co-Eruptive Modeling

We followed the approach of Lundgren et al. (2013) to model the co-eruptive deformation source. The modeling occurred in two steps: (1) we solve for a deflating magma reservoir beneath Taal using a CDM and a simple planar tensile dislocation (dike) model with uniform opening, and (2) after we fixed the parameters of the CDM and the rectangular dislocation, we developed a distributed opening model to further constrain the dike dimensions and opening using non-negative least squares method.

We used the Sentinel-1 InSAR ascending and descending cumulative LOS displacements and the descending range-offset maps for the co-eruptive modeling since these data provided better constraints on the initial activity of the dike. We downsampled each data set and calculated their covariances using a model-based quadtree approach (Lohman & Simons, 2005) given a dike model.

2.2.3. Post-Eruptive Modeling

We only modeled the deformation from January 15 to February 4, 2020 to track the further activity of the dike after the main eruptive phase using the Sentinel-1 ascending and descending cumulative LOS displacements. Following the steps in Section 2.2.2, we first fixed the geometries and locations of the magma reservoir and the dike either by: Case (1) using the co-eruptive MAP solutions (Table S1) or Case (2) independently estimating the dike parameters using a Markov chain Monte Carlo (MCMC) Bayesian inference method (Table S1). Afterward, we estimated the distributed opening of the dike as in step-2 of the co-eruptive modeling.
### 3. Crisis Response Story

During the peak of the eruption, monitoring instruments installed on Volcano Island were either destroyed or temporarily stopped operating (Sabillo, 2020). SAR data delivered with low-latency were essential in understanding the state of the volcano and in guiding PHIVOLCS to forecast possible eruption scenarios as the event progressed. On 12 January, we initiated a crisis response by producing interferograms using available datasets at that time: (1) ALOS-2 descending SAR datasets spanning one year (i.e., January 13, 2019–January 12, 2020) and (2) Sentinel-1 ascending SAR datasets covering 12 days prior to the eruption (i.e., December 30, 2019–January 11, 2020), to determine if there was any long- or short-term pre-eruptive deformation detected by InSAR. The ALOS-2 descending interferogram (Figure 2a) revealed roughly two fringes crossing Volcano Island (i.e., ~24 cm of LOS surface deformation) over one year, whereas the Sentinel-1 ascending interferogram (Figure 2b) showed null to very small short-term deformation. Regions outside the European territory are typically imaged by the Sentinel-1 satellite every 12 days. To monitor the on-going volcanic crisis with low-latency, the International Disaster Charter was triggered and, at the same time, we requested a 6-day repeat Sentinel-1 acquisition from the ESA, which operates Sentinel-1 for the Copernicus program of the European Commission. ESA granted the request until 22 February 2020.

On 15 January, after the descending SAR data became available in the Copernicus sci-hub repository (https://scihub.copernicus.eu/), we rapidly processed and directly communicated the first co-eruptive Sentinel-1 interferogram (Figure 3a) to PHIVOLCS through a messaging application—just in time for their science discussion and planning. The interferogram revealed dense fringes within Volcano Island and the Pansipit River Valley (PRV). However, due to the large deformation caused by the eruption, it was difficult to measure the ground displacements within these areas. The pixel-offset analysis (Figure 3c) enabled us to see through the regions that have significantly changed, however, with a lower accuracy compared to InSAR-derived phase measurements. We interpreted the signal from this first co-eruptive InSAR and pixel-offset maps as a result of the emplacement of a large dike and the deflation of a magma reservoir. This initial interpretation during the on-going crisis was further supported by the subsequent ascending interferogram and modeling (see Section 5). If a dike-fed eruption had occurred at that time, we would have inferred the dike breaking to the surface near the vicinity of the Pansipit River Valley, where the towns of Taal, Lemery, Agoncillo, and San Nicolas are directly affected. The information derived from low-latency InSAR data, coupled with the recorded seismicity during the eruption, enabled PHIVOLCS to monitor the volcanic state despite loss of several in situ instruments. We continued providing InSAR data and analyses to PHIVOLCS as the eruption progressed until it had completely waned.

### 4. The Pre-Eruptive State, Location, and Geometry of the Magma Storage Beneath Taal

The Sentinel-1 datasets stretch back from late 2014 until 1 day before Taal's 2020 eruption. For both the Sentinel-1 tracks (Figure 2c), the time-series revealed a negative trend in the measured LOS displacement from 2014 until late 2018 before the slope shifted to positive, leading up to the eruption day. This negative trend followed by a positive displacement trend has also been recorded through geodetic surveys conducted by PHIVOLCS (PHIVOLCS, 2019). The ALOS-2 time-series from both tracks started in late 2018 until a few hours before the eruption and thus only revealed the inflationary part of the deformation time-series.

From the inversion, we found a NE-SW striking ellipsoidal body located at ~5 km depth slightly east of Volcano Island (Figures S1–S2). Seismic analysis (Kumagai et al., 2014), magnetotelluric soundings (Zlotnicki et al., 2017), and geodetic studies of Taal volcano also showed a reservoir at roughly similar depth and location, although with a spherical-shaped chamber (Bartel et al., 2003; Galgana et al., 2014; Morales-Rivera et al., 2019; Zlotnicki et al., 2018). The estimated reservoir volume change was 0.045 ± 0.007 km³ over the year prior to the eruption (Figure S3a). Note that all our reported volume changes in this manuscript and their uncertainties are only based on simple assumptions and crustal rheology (i.e., elastic, homogenous half-space). Considering a non-elastic rheology will modify the volume estimation required to fit the surface deformation, however: (1) as shown by Holohan et al. (2017), taking into account the damage of the crust induced changes, which are significant only for a large depletion of the storage zone and for shallow reservoirs and (2) we lack constraints on the local rheology's deviation from an elastic one. The lack of constraint
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a and b) Pre-eruptive interferograms of Taal produced using (a) ALOS-2 L-band and (b) Sentinel-1 C-band satellite data set showing long- and short-term surface deformations, respectively. (c) Cumulative descending LOS displacement map from Sentinel-1 track 032, covering January 2, 2019–January 9, 2020. (d) Surface LOS displacement time-series measured at 14.0188° N, 121.0012° E (i.e., star symbol in Figure 2c) from both the ascending and descending tracks of ALOS-2 and Sentinel-1 starting in late 2014 and leading up to the January 12, 2020 Taal eruption (dashed red line). For each time-series, we referenced the date to their respective first acquisitions in January 2019 (ALOS-2 A138: 22 January 2019; ALOS-2 D027: January 13, 2019; S1 A142: January 4, 2019; S1 D032: January 2, 2019). Blue dots: displacement time-series from Sentinel-1 ascending track 142. Orange dots: displacement time-series from Sentinel-1 descending track 032. Orange circles: displacement time-series from Sentinel-1 descending track 032 generated using the NSBAS processing chain. Note that the NSBAS-derived result is only shown here for comparison and was not used during the inversion. Green triangles: displacement time-series from ALOS-2 track ascending track A138. Red triangles: displacement time-series from ALOS-2 descending track 027. The light blue and light orange lines correspond to the negative displacement trend for the Sentinel-1 derived measurements. LOS, line-of-sight.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
on the local rheology would significantly increase our estimated volume uncertainties. Figure S2 gives the data, model and residual maps using the MAP values illustrating a good fit between the model inversion results and data.

5. Magma Withdrawal and Dike Emplacement

From the co-eruptive InSAR phase (Figures 3a and 3b) and pixel-offset (Figures 3c and 3d) displacement maps, we derived the horizontal (Figure 3e) and vertical displacements (Figures 3f and 3g). Two distinct signals emerged: (1) Volcano Island and the north and east portions of Taal caldera experienced deflation of ~4 m, and (2) the SW region from Volcano Island toward Balayan Bay inflated (~1 m) and pulled apart (~2 m), with the center of the NE-SW trending rift located near the Pansipit river. In volcano geodesy studies, these signals suggest magma withdrawal from a reservoir and magma emplacement through dike intrusion (Neal et al., 2019; Shreve et al., 2019; Sigmundsson et al., 2015; Sturkell et al., 2006; Wright et al., 2006, & references therein).

From the co-eruptive modeling, we find that the magma reservoir at ~5-6 km depth below Volcano Island lost an estimated volume of -0.531 ± 0.004 km³ (Figures S3b and S4), which is an order of magnitude larger than the accrued volume change of the reservoir the prior year. This indicates that only a small percentage of the magma withdrawn from the reservoir came from the short-term accumulation period. The model also revealed that a 21 × 8 km, near-vertical, NE-striking dike was intruded below the surface from Volcano Island extending southwestward toward Balayan Bay (Figure S5). NE-SW trending ground fissures that emerged in the municipalities of Agoncillo, Lemery, San Nicolas, Taal, and Talisay, matched the results of our dike model (PHIVOLCSd, 2020). We estimated that the volume of the dike is 0.643 ± 0.001 km³ (Figure S3c) assuming a uniform opening, which has the same order of magnitude as the volume loss of the magma reservoir. The volume of the dike is only 1.2× greater than the volume loss of the magma reservoir, which can imply that the volume difference may be due to additional deep source input or the effect of magma compressibility as a result of gas exsolution (Rivalta & Segall, 2008). Fixing the magmatic source to a deeper depth (i.e., >6 km) might result in equal volumes, however there is currently no independent data set to support this claim. During the on-going eruption response, we estimated that if half of the dike's volume had erupted explosively, it might have yielded a Volcano Explosivity Index (VEI) 4 eruption (Newhall & Self, 1982).

The co-eruptive distributed opening model in Figures 3h–3i indicated strong opening below Volcano Island that extended southwestward toward the Pansipit River Valley at <10 km depth. Although the strongest opening is found beneath Taal Lake, it is unclear if the dike reached the bottom of the lake, which has a maximum depth of 198 m (Castillo & Gonzales, 1976). Notice that the co-eruptive model also shows opening at around 15–20 km depth indicating that there may be a possible deep source feeding the shallow reservoir, although this can also represent the inability of the shallower portions of the model to fit far-field LOS displacements that may be dominated by spatially long wavelength sources of noise in the InSAR data. The estimated volume of the dike from the distributed opening model is 0.689 km³, which has the same order of magnitude as the rectangular dislocation model. Figure S7 shows the co-eruptive data, model, and residual maps using the MAP values, illustrating a good fit between the model and the data.

Subsequent to the main eruptive phase, persistent surface deformation was still observed at Taal (Figures 4a–4d; Figures S8–S12). Volcano Island continued to deflate of up to 0.08 m until around 04 February (Figures 4a and 4c) before inflationary signals are detected (Figures 4b and 4d; Figures S8 and S10). Ground displacements are sustained around the Pansipit River Valley toward Balayan Bay, with most of the significant deformation occurring almost simultaneous to the deflation of Volcano Island. We observed that the

**Figure 3.** The co-eruptive datasets and model. (a and b) The InSAR phase, re-wrapped at 5 cm per color cycle, and (c and d) range-offset displacement maps from Sentinel-1 descending and ascending tracks showing how much the ground moved as a result of the Taal eruption and dike emplacement. The green rectangles in Figures 3a and 3b is the trace of the dike from the model inversion. (e and f) Derived horizontal and vertical displacements from InSAR phase and range-offset maps. Note that areas covered by fringe aliasing in the InSAR maps are masked and are replaced by range-offset data. (g) Magnified view of the vertical displacement map to highlight the local subsidence within the island and the location of emplaced volcanic deposits (encircled in red). (h) Map view and (i) 3D view of the CDM and distributed opening model. PRV: Pansipit River Valley; VI: Volcano Island. (j) Solution root-mean-square misfit versus roughness as a function of the smoothing value used for the dike. The value in red text is the smoothing factor that we used for the co-eruptive model.
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Post-eruptive data. Horizontal (a and b) and vertical displacement (c and d) maps derived using Sentinel-1 descending and ascending datasets covering two periods: (a and c) 15 January to 04 February and (b and d) 02 February to 27 June 2020. (a and b) Red means that the ground moved eastward, blue signifies that the ground moved westward. (c and d) Red means that the ground moved upward and blue represents downward movement.</figcaption></figure>

reservoir continued to lose volume (i.e., -0.00775 ± 0.0001 km³, Figure S3d) while the dike continued to grow (i.e., 0.0557 ± 0.0004 km³, Figure S3e) until 4 February.

The post-main eruptive phase distributed opening model between 15 January and 4 February shows weak opening (i.e., up to 41 cm) centered mostly above 10 km-depth beneath Taal lake for both the cases described in Section 2.2.3 (Figures S13 and S17). No further dike propagation is also evidenced, although the resolution of the dike tip position is limited by the absence of information on the displacement field under the sea. The model in Figure S13 follows the conditions in Case 1 and can only partially explain the post-eruptive data (Figure S14). We reduced the residuals (Figure S18) only when all the parameters of the dike were re-estimated (Figure S17) prior to Step 2, implying that the dike's orientation may have changed as it propagated (i.e., shift to a more vertically dipping dike, Figure S15-S18). There are also transient processes (i.e., local faulting events in the graben area and the emergence of fissures in the northern part of the caldera [PHIVOLCSd, 2020]) that were not accounted for in the post-eruptive model.
Between 2 February and 27 June, InSAR time-series of Taal reveal that (1) the magma reservoir beneath Volcano Island is in recovery as evidenced by the inflation detected around the island, and (2) no significant surface deformation was detected in the Pansipit River Valley and in the nearby areas (Figures S8-S12). The re-inflation and thus, the refilling of the magma reservoir is consistent with our co-eruptive model showing possible evidence of a shallow reservoir fed by a deeper source.

## 6. Discussion

When a volcanic eruption occurs, we ask the obvious however difficult fundamental questions. *What drives the eruption? Why and when did the shallow magmatic storage become unstable and fail?* And for Taal's specific case, *why did the magma favor lateral propagation and remain stalled at depth?*

Magma produced in the mantle is thought to rise in the form of dike or melt-pocket injections through cracks in the crust driven by its buoyancy relative to the surrounding crust. It can stall at depth once it loses buoyancy or cools and sometimes when it encounters a barrier. Magma may accumulate at shallower zones before erupting to the surface, form magma reservoirs through series of intrusions, and potentially displaces and interacts with older magmas (NAP, 2017).

In the upper 10 km of the crust, the controls on explosivity and eruption include a wide range of intrinsic magma properties (i.e., magma ascent rate, volatile content, magma rheology/viscosity, magma chamber overpressure, magma temperature, magma porosity and permeability, magma depth/pressure, crystal and bubble content) and extrinsic parameters (i.e., stress field changes, decompression rate, crustal properties, conduit and vent geometry, external assimilation of surface water/carbonates) that have competing effects and potentially contradictory outcomes (Cassidy et al., 2018). Nearer the surface, eruption explosivity becomes highly dependent on the initial gas content and melt and magma viscosity (Newhall, 2015). For this reason, high-gas content and high-viscosity magmas often result in explosive eruptions.

How eruptions are triggered is a complex interaction of magma influx into the reservoir (recharge), the size of the reservoir, and the gas content of the magma (NAP, 2017; Degruyter & Huber, 2014, & references therein). Recently, a model whereby passive degassing causes decompression (deflation) followed by rapid magma ascent and pressurization (inflation) of the shallow reservoir has been proposed as a possible mechanism for eruption initiation at volcanoes with hydraulically connected plumbing systems (Girona et al., 2014, 2015). In this model, degassing creates a pressure imbalance that promotes magma ascent from the deep to shallow magma reservoirs. Time-series of deflation-inflation events over multi-year time scales have been observed at numerous systems both with (e.g., Copahue, Lundgren et al., 2017; Velez et al., 2016, 2011) and without eruptions (e.g., Campi Flegrei, Troise et al., 2019; Long Valley, Montgomery-Brown et al., 2015; Domuyo, Lundgren et al., 2020). In the case of Copahue volcano, the linear deflation from 2004 to 2008 was attributed to leakage of magmatic fluids through a self-sealing carapice overlying a multiphase brine surrounding a shallow magma body, whereas the inflation from 2011 to 2016 was found related to the transfer of magma in a two-source system.

As the magma reservoir pressurizes, it may rupture and initiate the propagation of a dike once a critical overpressure is reached. The balance between vertical and lateral propagation is a key question central to many magmatic systems, occurring not only in tectonic spreading centers, such as Krafla, Afar, and Bárdarbunga (Sigmundsson et al., 2015; Sturkell et al., 2006; Wright et al., 2006), but also at ocean islands with localized rift zones such as Kilauea, Ambrym, and Piton de la Fournaise (Froger et al., 2015; Neal et al., 2019; Shreve et al., 2019) and sometimes even at restless calderas such as Campi Flegrei (Di Vito et al., 2016). In these systems, a large volume of magma travels laterally, whether the eruptive sequence starts with an eruption in the central area, as observed at Piton de la Fournaise in 2007 (Froger et al., 2015) and Ambrym in 2018 (Shreve et al., 2019), or without eruption in the central part but deflation due to the lateral dike propagation and distal eruption as found for Bárdarbunga-Holuhraun, 2015 (Sigmundsson et al., 2015) and Kilauea, 2018 (Neal et al., 2019). However, there are also cases where lateral dike propagation does not feed an eruption (Sturkell et al., 2006; Wright et al., 2006), as in the case of Taal's 2020 event.

The direction of magma propagation is mainly controlled by the local stress field and the magma driving pressure due to buoyancy and bottom magma influx, (Pinel et al., 2017) as well as local rheology (Urbani
et al., 2017, 2018). Indeed, positive buoyancy will always favor magma ascent and in all models of lateral propagation (Buck et al., 2006; Einarsson et al., 1980; Grandin et al., 2012; Heimisson et al., 2015; Pinel & Jaupart, 2004; Townsend et al., 2017; Urbani et al., 2017, & references therein), buoyancy should act against the direction of propagation being controlled by the local stress field and/or the topography.

For Taal's case, its ~5–6 km deep shallow reservoir became unstable and ruptured due to a small magma input during the short-term accumulation period. Taal has a history of frequent seismovolcanic crises and/ or surface activity which are interpreted as repeated magma recharge of its shallow reservoir (Zlotnicki et al., 2018). The sudden magma injection in 2019 might have been due to the pressure imbalance induced by the depressurization of the shallow storage zone, resulting from long-term diffuse degassing, as supported by: (1) the sequence of long-term deflation and short-term inflation detected by InSAR, and (2) the decrease in CO₂ flux measurements from 2015 to 2017 before rapidly rising from 700 T/day (background value) to 2400 T/day between February 2019 and the eruption day (Bernard et al., 2020).

As the magma reservoir accumulated pressure, it also interacted with the overlying hydrothermal system, resulting to the initial phreatic eruption (Figure 1c) which potentially led to the failure of the shallow reservoir. The timing of events can be further constrained with the availability of seismic, GNSS, and petrological datasets. Nevertheless, current knowledge of Taal's plumbing system from geochemical and geophysical studies (Alanis et al., 2013; Maussen et al., 2018; Yamaya et al., 2013; Zlotnicki et al., 2018), and CO₂ flux and temperature measurements from Main Crater Lake (Bernard et al., 2020), indicate that high-flux of gases and volcanic fluids are transferred from the shallow magma reservoir to the hydrothermal system (Zlotnicki et al., 2018). This eventually leads to three possible scenarios, the hydrothermal system may: (1) have been heated, sealed and pressurized, (2) have accumulated magmatic gas resulting in a reduced boiling point and overpressurization of the system (Bernard et al., 2020), or 3) have been destabilized due to the injection of magma, either as small batches or through interaction with a large dike. It is worth mentioning that the volume of erupted tephra (i.e., ~0.032 km³, VEI = 3) from 12–13 January 2020 (Martinez-Villegas, 2020), is close to the estimated: (1) volume of products that can be initially involved if a phreatic eruption occurs at Taal due to the mechanically weak and mineralized materials in the northern part of Volcano Island (Zlotnicki et al., 2018) and (2) pre-eruptive volume increase of the magma reservoir from 2019.

The phreatic-phreatomagmatic eruption that occurred in Volcano Island was driven by an initial gas-rich phase (Bernard et al., 2020) which further decreased the remaining gas in the shallow storage zone. The orientation of the lateral dike is also consistent with the pre-existing fault structures along the NE-SW-striking Macolod Corridor extensional zone, and appear to have played a significant role during the diking event by promoting the opening in the direction of the minimum compressive stress and the propagation along the maximum compressive stress (Anderson, 1951; Cotterell & Rice, 1980; Delaney et al., 1986). The fact that Taal's topography is very low could have also favored the deflection of the magma (Corbi et al., 2015; Gaete et al., 2019). Both the deflection of the magma path toward the edge of the caldera and the predominance of a lateral transport compared to the vertical one indicate that the magma buoyancy was reduced, probably due to the efficient gas loss at pre-eruptive and initial eruptive phases. This low gas content is also consistent with the fact that the dike volume is only slightly larger than the volume loss from the reservoir.

Although our model showed possible rupture in the lake between the Pansipit River Valley and Volcano Island, there is no reported evidence of underwater lava flow emplacement and our datasets are restricted inland such that we do not have enough constraints beneath the water. Post-eruption bathymetric studies and geochemical analyses of Taal lake will provide more information and constraints. Our models indicate that the dike migrated at deeper depths toward the Pansipit River Valley and the post-eruptive dike growth was limited to a widening beneath Taal lake. Such observations are in agreement with dike propagation models controlled by local topography, thus inhibiting vertical propagation (Maccaferri et al., 2016).

Combining all these pieces of evidence (i.e., gas loss during pre-eruptive and initial eruptive phase, presence of pre-existing weaknesses, reduced magma buoyancy), we conclude that the Taal 2020 volcanic crisis did not result in an eruption, despite its history and potential of creating large and destructive events, simply because there was insufficient pressure to drive the magma upward and hence it remained stalled at depth. Eruptions in the Main Crater and ground fissuring related to dike intrusions or fault adjustments were reported in the past (Delos Reyes et al., 2018). Although the 2020 crisis was not highly destructive, the past
eruptions closest to this event were probably those in 1749 and 1911, which were categorized as VEI 3–4 events (Delos Reyes et al., 2018). Both eruptions were characterized by violent phreatomagmatic or plinian eruptions in the Main Crater while NE-SW-striking ground fissures appeared (1) at the northern part of the volcano reaching as far as Calamba (1749) or (2) from Lemery to Calamba (1911) (Delos Reyes et al., 2018).

## 7. Lessons Learned

Taal is a closely monitored volcano, consists of multi-parametric network of ground instruments. This is the first time that comprehensive and well-documented InSAR-derived displacement maps, time-series analyses, and deformation modeling of Taal's pre-, co-, and post-eruptive state are presented.

On one hand, the pre-eruptive InSAR time-series revealed the inflation during the year preceding the January 2020 eruption. Although this inflation has also been detected by the geodetic network of Taal, InSAR provided the complementary spatial information. We obtained the geometry, location, and source strength of the magma reservoir beneath Volcano Island which are helpful information when characterizing, in our case, the minimum potential magnitude of an impending eruption. On the other hand, the co-eruptive interferogram showed the emplacement of the large dike allowing us to constrain the possible location of a dike-fed eruption had it occurred at that time and to immediately estimate the maximum VEI, which are important information for crisis management.

Although satellite data will never replace in situ observations, they can provide complementary and significant information with unprecedented spatial coverage. Indeed, multi-parametric datasets that are comprised of ground- and satellite-based measurements should be used in monitoring volcanoes to get a full perspective of the volcanic state and to obtain insights into poorly known volcanic processes.

The launch of the Sentinel-1 satellites in 2014 and 2016 has greatly improved InSAR science and applications, particularly, in the monitoring of crustal deformation, allowing broad access to free and low-latency SAR data. Taal's recent crisis greatly benefited from the freely available Sentinel-1 data and the continuous and real-time delivery of processed data, analyses, and models that were crucial to the observatory when most ground instruments were no longer available. During the eruption, the biggest challenge that we encountered was related to how accurate and correct were the low-latency preliminary analyses and models that we provided to PHIVOLCS, given the lack of prior InSAR baseline observations for Taal. Fortunately, although the results that we report in this manuscript required extensive reprocessing and utilization of larger datasets and tools, our analyses remained consistent with those delivered to the observatory during the height of the crisis and the overall behavior of the volcano also corroborated our results. Despite the unavailability of other co-eruptive in situ data, actual surface deformation observed in the field such as the emergence of ground fissures and the occurrence of subsidence and uplift in the coastal lake regions based on available field photographs confirmed our analyses and interpretations. We, therefore, highlight the need and encourage volcano observatories globally to include remote sensing data like InSAR as a major component in their regular monitoring of active volcanoes, particularly now that (1) the data and processing tools are becoming more open and freely available, and (2) advanced satellite system acquisitions have become more suitable for highly vegetated regions.

## Data Availability Statement

All the ARIA-produced interferograms used in this study are accessible via the ASF DAAC vertex page (https://search.asf.alaska.edu). All the ALOS-2 interferograms are available on Zenodo at https://doi.org/10.5281/zenodo.4574537.

## References

* Alanis, P. K., Yamaya, Y., Takeuchi, A., Sasai, Y., Okada, Y., & Nagao, T. (2013). A large hydrothermal reservoir beneath taal volcano (philippines) revealed by magnetotelluric observations and its implications to the volcanic activity. *Proceedings of the Japan Academy, Series B*, 89(8), 383–389.

* Anderson, E. M. (1951). *The dynamics of faulting and dyke formation with applications to Britain*. Oliver and Boyd.

* Aurelio, M. A (2000). Shear partitioning in the Philippines: Constraints from Philippine Fault and global positioning system data. *Island Arc*, 9(4), 584–597. https://doi.org/10.1046/j.1440-1738.2000.00304.x

## Acknowledgments

We thank Hook Hua, Lan Dang, and David Bekaert for activating and managing the ARIA Project Sentinel-1A/B SAR interferogram generation based on European Commission Copernicus Sentinel-1 satellite SAR data products courtesy of the European Space Agency. We thank the Japan Aerospace Exploration Agency (JAXA) for providing ALOS-2 PALSAR-2 SAR L1.1 (SLC) data through RA6 Proposal P3024002 (PI Lundgren) used in the analyses; original L1.1 products ©JAXA 2019–2020. We are also grateful to Chris Newhall and Eric Fielding for the helpful discussions during the peak of the eruption. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (grant 281945.02.47.04.51). AlTar software development and modeling supported through the JPL R&TD program (project number 01STCR- R.18.245.034). V. Pinel acknowledges support through the MagmaPropagator project (contract ANR-18-CE92-0037-01). We thank Maurizio Battaglia and our two anonymous reviewers for their constructive suggestions that improved our manuscript.
Bartel, B. A., Hamburger, M. W., Meertens, C. M., Lowry, A. R., & Corpuz, E. (2003). Dynamics of active magmatic and hydrothermal systems at Taal Volcano, Philippines, from continuous GPS measurements. *Journal of Geophysical Research*, **108**(B10). https://doi.org/10.1029/2002jb002194

Bekaert, D., Karim, M., Justin, L., Hua, H., Agram, P., Owen, S., & others (2019). Development and dissemination of standardized geodetic products by the advanced rapid imaging and analysis (ARIA) center for natural hazards. In *Proceedings of the International Union of Geodesy and Geophysics (IUGG)* (pp. 8–18). Montreal, Canada.

Bernard, A., Villacorte, E., Maussen, K., Caudron, C., Robic, J., Maximo, R., & Solidum, R. (2020). Carbon dioxide in Taal volcanic lake: A simple gasometer for volcano monitoring. *Geophysical Research Letters*, **47**(24), e2020GL090884. https://doi.org/10.1029/2020GL090884

Buck, W. R., Einarsson, P., & Brandsdóttir, B. (2006). Tectonic stress and magma chamber size as controls on dike propagation: Constraints from the 1975-1984 Krafla rifting episode. *Journal of Geophysical Research*, **111**(B12). https://doi.org/10.1029/2005jb003879

Buzzanga, B., Bekaert, D. P., Hamlington, B. D., & Sangha, S. S. (2020). Toward sustained monitoring of subsidence at the coast using InSAR and GPS: An application in hampton roads, Virginia. *Geophysical Research Letters*, **47**(18), e2020GL090013. https://doi.org/10.1029/2020gl090013

Cassidy, M., Manga, M., Cashman, K., & Bachmann, O. (2018). Controls on explosive-effusive volcanic eruption styles. *Nature Communications*, **9**(1), 1–16. https://doi.org/10.1038/s41467-018-05293-3

Castillo, B., & Gonzales, C. (1976). Hydrology of Taal Lake. *Fisheries Research Journal of the Philippines*, **1**(2), 62–75.

Corbi, F., Rivalta, E., Pinel, V., Maccaferri, F., Bagnardi, M., & Acocella, V. (2015). How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes. *Earth and Planetary Science Letters*, **431**, 287–293. https://doi.org/10.1016/j.epsl.2015.09.028

Cotterell, B., & Rice, J. R. (1980). Slightly curved or kinked cracks. *International Journal of Fracture*, **16**(2), 155–169. https://doi.org/10.1007/bf00012619

Degruyter, W., & Huber, C. (2014). A model for eruption frequency of upper crustal silicic magma chambers. *Earth and Planetary Science Letters*, **403**, 117–130. https://doi.org/10.1016/j.epsl.2014.06.047

Delaney, P. T., Pollard, D. D., Ziony, J. I., & McKee, E. H. (1986). Field relations between dikes and joints: Emplacement processes and paleostress analysis. *Journal of Geophysical Research*, **91**(B5), 4920–4938. https://doi.org/10.1029/jb091ib05p04920

Delos Reyes, P. J., Bornas, M. A. V., Dominey-Howes, D., Pidlaoan, A. C., Magill, C. R., & Solidum, Jr, R. U. (2018). A synthesis and review of historical eruptions at Taal Volcano, Southern Luzon, Philippines. *Earth-Science Reviews*, **177**, 565–588. https://doi.org/10.1016/j.earscirev.2017.11.014

De Zan, F., Ansari, H., Gomba, G., Brcic, R., & Parizzi, A. (2019). *The impact of closure phases on InSAR processing*. In Agu fall meeting abstracts (Vol. 2019). G21A-05.

Di Vito, M. A., Acocella, V., Aiello, G., Barra, D., Battaglia, M., Carandente, A., & others (2016). Magma transfer at Campi Flegrei Caldera (Italy) before the 1538 AD eruption. *Scientific Reports*, **6**(1), 1–9. https://doi.org/10.1038/srep32245

Doin, M.-P., Guillaaso, S., Jolivet, R., Lasserre, C., Lodge, F., Ducret, G., & Grandin, R. (2011). Presentation of the small baseline NSBAS processing chain on a case example: The Etna deformation monitoring from 2003 to 2010 using Envisat data. In *Proceedings of the fringe symposium* (pp. 3434–3437).

Einarsson, P., & Brandsdóttir, B. (1980). Seismological evidence for lateral magma intrusion during the July 1978 deflation of the Krafla volcano in NE-Iceland. *Journal of Geophysics*, **47**(1), 160–165.

Fikos, I., Vargemezis, G., Zlotnicki, J., Puertollano, J. R., Alanis, P. B., Pigtain, R. C., et al. (2012). Electrical resistivity tomography study of Taal volcano hydrothermal system, Philippines. *Bulletin of Volcanology*, **74**, 1821–1831. https://doi.org/10.1007/s00445-012-0638-5

Friger, J.-L., Famin, V., Cayol, V., Augier, A., Michon, L., & Lénat, J.-F. (2015). Time-dependent displacements during and after the April 2007 eruption of Piton de la Fournaise, revealed by interferometric data. *Journal of Volcanology and Geothermal Research*, **296**, 55–68. https://doi.org/10.1016/j.jvolgeores.2015.02.014

Gaete, A., Kavanagh, J. L., Rivalta, E., Hilmi Hazim, S., Walter, T. R., & Dennis, D. J. C. (2019). The impact of unloading stresses on post-caldera magma intrusions. *Earth and Planetary Science Letters*, **508**, 109–121. https://doi.org/10.1016/j.epsl.2018.12.016

Galgana, G. A., Newman, A. V., Hamburger, M. W., & Solidum, R. U. (2014). Geodetic observations and modeling of time-varying deformation at Taal Volcano, Philippines. *Journal of Volcanology and Geothermal Research*, **271**, 11–23. https://doi.org/10.1016/j.jvolgeores.2013.11.005

Girona, T., Costa, F., Newhall, C., & Taisne, B. (2014). On depressurization of volcanic magma reservoirs by passive degassing. *Journal of Geophysical Research: Solid Earth*, **119**(12), 8667–8687. https://doi.org/10.1002/2014jb011368

Girona, T., Costa, F., & Schubert, G. (2015). Degassing during quiescence as a trigger of magma ascent and volcanic eruptions. *Scientific Reports*, **5**, 18212. https://doi.org/10.1038/srep18212

Grandin, R. (2015). *Interferometric processing of SLC Sentinel-1 TOPS data*.

Grandin, R., Socquet, A., Doubre, C., Jacques, E., & C.P. King, G. (2012). Elastic thickness control of lateral dyke intrusion at mid-ocean ridges. *Earth and Planetary Science Letters*, **319–320**, 83–95. https://doi.org/10.1016/j.epsl.2011.12.011

Heimisson, E. R., Hooper, A., & Sigmundsson, F. (2015). Forecasting the path of a laterally propagating dike. *Journal of Geophysical Research: Solid Earth*, **120**(12), 8774–8792. https://doi.org/10.1002/2015jb012402

Holohan, E. P., Sudhaus, H., Walter, T. R., Schöpfer, M. P., & Walsh, J. J. (2017). Effects of host-rock fracturing on elastic-deformation source models of volcano deflation. *Scientific Reports*, **7**(1), 1–12. https://doi.org/10.1038/s41598-017-10009-6

Kumagai, H., Lacson, R., Jr, Maeda, Y., Figueroa, M. S., & Yamashina, T. (2014). Shallow S wave attenuation and actively degassing magma beneath Taal Volcano, Philippines. *Geophysical Research Letters*, **41**(19), 6681–6688. https://doi.org/10.1002/2014gl061193

Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. *Geochemistry, Geophysics, Geosystems*, **6**(1). https://doi.org/10.1029/2004gc000841

Lundgren, P., Girona, T., Bato, M. G., Realmuto, V. J., Samsonov, S., Cardona, C., & Aivazis, M. (2020). The dynamics of large silicic systems from satellite remote sensing observations: The intriguing case of Domuyo volcano, Argentina. *Scientific Reports*, **10**(1), 1–15. https://doi.org/10.1038/s41598-020-67982-8

Lundgren, P., Nikkhoo, M., Samsonov, S. V., Milillo, P., Gil-Cruz, F., & Lazo, J. (2017). Source model for the copahue volcano magma plumbing system constrained by InSAR surface deformation observations. *Journal of Geophysical Research: Solid Earth*, **122**(7), 5729–5747. https://doi.org/10.1002/2017jb014368

Lundgren, P., Poland, M., Miklius, A., Orr, T., Yun, S.-H., Fielding, E., et al. (2013). Evolution of dike opening during the March 2011 Kamoamoa fissure eruption, Kilauea Volcano, Hawaii. *Journal of Geophysical Research: Solid Earth*, **118**(3), 897–914. https://doi.org/10.1002/jgrb.50108

Maccaferri, F., Rivalta, E., Passarelli, L., & Aoki, Y. (2016). On the mechanisms governing dike arrest: Insight from the 2000 Miyakejima dike injection. *Earth and Planetary Science Letters*, **434**, 64–74. https://doi.org/10.1016/j.epsl.2015.11.024
Martinez-Villegas, M. (2020). *Volcano risk communication in the Philippines: Challenges, lessons and opportunities Taal volcano 2020*. Retrieved from https://drive.google.com/file/d/1sUrKHfObZKDBA8y1sDINr9\text[\_\]CcqPORSWB

Maussen, K., Villacorte, E., Rebadulla, R. R., Maximo, R. P., Debaille, V., Bornas, M. A., & Bernard, A. (2018). Geochemical characterization of Taal volcano-hydrothermal system and temporal evolution during continued phases of unrest (1991-2017). *Journal of Volcanology and Geothermal Research*, **352**, 38–54. https://doi.org/10.1016/j.jvolgeores.2018.01.007

Minson, S. E., Simons, M., & Beck, J. L. (2013). Bayesian inversion for finite fault earthquake source models I-theory and algorithm. *Geophysical Journal International*, **194**(3), 1701–1726. https://doi.org/10.1093/gji/ggt180

Montgomery-Brown, E. K., Wicks, C. W., Cervelli, P. F., Langbein, J. O., Svarc, J. L., Shelly, D. R., et al. (2015). Renewed inflation of long valley caldera, California (2011 to 2014). *Geophysical Research Letters*, **42**(13), 5250–5257. https://doi.org/10.1002/2015gl064338

Morales Rivera, A. M., Amelung, F., Albino, F., & Gregg, P. M. (2019). Impact of crustal rheology on temperature-dependent viscoelastic models of volcano deformation: Application to Taal volcano, Philippines. *Journal of Geophysical Research: Solid Earth*, **124**(1), 978–994. https://doi.org/10.1029/2018jb016054

NAP. (2017). *Volcanic eruptions and their repose, unrest, precursors, and timing*. National Academies Press.

NDRRMC. (2000). *Situational report No. 74 re Taal volcano eruption*. Retrieved from https://reliefweb.int/sites/reliefweb.int/files/resources/SitRep_No_74_re_Taal_Volcano_Eruption_issued_on_13_Feb_2020_6AM.pdf

Neal, C. A., Brantley, S. R., Antolik, L., Babb, J. L., Burgess, M., Calles, K., et al. (2019). The 2018 rift eruption and summit collapse of Kilauea Volcano. *Science*, **363**(6425), 367–374. https://doi.org/10.1126/science.aav7046

Newhall, C. (2015). 4.14—Volcanology 101. In G. Schubert (Ed.), *Treatise on geophysics* (2nd ed., pp. 355–387). Oxford, UK: Elsevier. https://doi.org/10.1016/B978-0-444-53802-4.00080-4

Newhall, C. G., & Self, S. (1982). The volcanic explosivity index (VEI) an estimate of explosive magnitude for historical volcanism. *Journal of Geophysical Research*, **87**(C2), 1231–1238. https://doi.org/10.1029/jc087ic02p01231

Nikkhoo, M., Walter, T. R., Lundgren, P. R., & Prats-Iraola, P. (2017). Compound dislocation models (CDMs) for volcano deformation analyses. *Geophysical Journal International*, **208**(2), 877–894.

PHIVOLCS. (2019). *Taal volcano bulletin 28 March 2019 9:15 am*. Retrieved from https://www.phivolcs.dost.gov.ph/index.php/taal-volcano-bulletin-menu/8117-taal-volcano-bulletin-28-march-2019-9-15-am-2

PHIVOLCSa. (2020). *Taal volcano bulletin 12 January 2020 07:30 pm*. Retrieved from https://www.phivolcs.dost.gov.ph/index.php/taal-volcano-bulletin-menu/9620-taal-volcano-bulletin-12-january-2020-07-30-pm

PHIVOLCSd. (2020). *General location of fissures related to the January 2020 Taal volcano eruptive activity*. Retrieved from https://www.phivolcs.dost.gov.ph/index.php/news/9662-general-location-of-fissures-related-to-the-january-2020-taal-volcano-eruptive-activity

Pinel, V., Carrara, A., Maccaferri, F., Rivalta, E., & Corbi, F. (2017). A two-step model for dynamical dike propagation in two dimensions: Application to the July 2001 Etna eruption. *Journal of Geophysical Research: Solid Earth*, **122**(2), 1107–1125. https://doi.org/10.1002/2016jb013630

Pinel, V., & Jaupart, C. (2004). Magma storage and horizontal dyke injection beneath a volcanic edifice. *Earth and Planetary Science Letters*, **221**(1–4), 245–262. https://doi.org/10.1016/s0012-821x(04)00076-7

Pubellier, M., Garcia, F., Loevenbruck, A., & Chorowicz, J. (2000). Recent deformation at the junction between the North Luzon block and the Central Philippines from ERS-1 Images. *Island Arc*, **9**(4), 598–610. https://doi.org/10.1046/j.1440-1738.2000.00305.x

Rimando, R. E., & Knuepfer, P. L. (2006). Neotectonics of the Marikina valley fault system (MVFS) and tectonic framework of structures in northern and central Luzon, Philippines. *Tectonophysics*, **415**(1–4), 17–38. https://doi.org/10.1016/j.tecto.2005.11.009

Rivalta, E., & Segall, P. (2008). Magma compressibility and the missing source for some dike intrusions. *Geophysical Research Letters*, **35**(4). https://doi.org/10.1029/2007gl032521

Rosen, P. A., Gurrola, E. M., Agram, P., Cohen, J., Lavalle, M., Riel, B., & Buckley, S. M. (2018). The InSAR scientific computing environment 3.0: A flexible framework for NISAR operational and user-led science processing. In *IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium* (pp. 4897–4900). Valencia: IEEE.

Sabillo, K. (2020). *Phivolcs cleans up ash-covered instruments, sets up new stations near Taal volcano*. Retrieved from https://news abs-cbn. com/news/01/15/20/phivolcs-cleans-up-ash-covered-instruments-sets-up-new-stations-near-taal-volcano

Shreve, T., Grandin, R., Boichu, M., Garaebiti, E., Moussallam, Y., Ballu, V., & others (2019). From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): The influence of regional tectonics. *Scientific Reports*, **9**(1), 1–13. https://doi.org/10.1038/s41598-019-55141-7

Sigmundsson, F., Hooper, A., Hreinsdóttir, S., Vogfjörð, K. S., Ófeigsson, B. G., Heimisson, E. R., et al. (2015). Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland. *Nature*, **517**(7533), 191–195. https://doi.org/10.1038/nature14111

Sturkell, E., Einarsson, P., Sigmundsson, F., Geirsson, H., Olafsson, H., Pedersen, R., & Stefánsson, R. (2006). Volcano geodesy and magma dynamics in Iceland. *Journal of Volcanology and Geothermal Research*, **150**(1–3), 14–34. https://doi.org/10.1016/j.jvolgeores.2005.07.010

Townsend, M. R., Pollard, D. D., & Smith, R. P. (2017). Mechanical models for dikes: A third school of thought. *Tectonophysics*, **703–704**, 98–118. https://doi.org/10.1016/j.tecto.2017.03.008

Troise, C., De Natale, G., Schiavone, R., Somma, R., & Moretti, R. (2019). The Campi Flegrei caldera unrest: Discriminating magma intrusions from hydrothermal effects and implications for possible evolution. *Earth-Science Reviews*, **188**, 108–122. https://doi.org/10.1016/j.earscirev.2018.11.007

Urbani, S., Acocella, V., & Rivalta, E. (2018). What drives the lateral versus vertical propagation of dikes? Insights from analogue models. *Journal of Geophysical Research: Solid Earth*, **123**(5), 3680–3697. https://doi.org/10.1029/2017jb015376

Urbani, S., Acocella, V., Rivalta, E., & Corbi, F. (2017). Propagation and arrest of dikes under topography: Models applied to the 2014 Bardarbunga (Iceland) rifting event. *Geophysical Research Letters*, **44**(13), 6692–6701. https://doi.org/10.1002/2017gl073130

Velez, M., Euillades, P., Blanco, M., & Euillades, L. (2016). Ground deformation between 2002 and 2013 from InSAR observations. In *Copahue volcano* (pp. 175–198). Springer.

Velez, M., Euillades, P., Caselli, A., Blanco, M., & Díaz, J. (2011). Deformation of copahue volcano: Inversion of InSAR data using a genetic algorithm. *Journal of Volcanology and Geothermal Research*, **202**(1–2), 117–126. https://doi.org/10.1016/j.jvolgeores.2011.01.012

Wolfe, J. A., & Self, S. (1983). Structural lineaments and Neogene volcanism in southwestern Luzon. *GMS*, **27**, 157–172. https://doi.org/10.1029/gm027p0157

Wright, T. J., Ebinger, C., Biggs, J., Ayele, A., Yirgu, G., Keir, D., & Stork, A. (2006). Magma-maintained rift segmentation at continental rupture in the 2005 Afar dyking episode. *Nature*, **442**(7100), 291–294. https://doi.org/10.1038/nature04978

Wright, T. J., Parsons, B. E., & Lu, Z. (2004). Toward mapping surface deformation in three dimensions using InSAR. *Geophysical Research Letters*, **31**(1). https://doi.org/10.1029/2003gl018827
Yamaya, Y., Alanis, P., Takeuchi, A., Cordon, J., Mogi, T., Hashimoto, T., & Nagao, T. (2013). A large hydrothermal reservoir beneath Taal volcano (Philippines) revealed by magnetotelluric resistivity survey: 2d resistivity modeling. *Bulletin of Volcanology*, 75(7), 1–13. https://doi.org/10.1007/s00445-013-0729-y

Yunjun, Z., Fattahi, H., & Amelung, F. (2019). Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction. *Computers & Geosciences*, 133, 104331. https://doi.org/10.1016/j.cageo.2019.104331

Zlotnicki, J., Sasai, Y., Johnston, M., Fauquet, F., Villacorte, E., & Cordon, J. (2018). The 2010 seismovolcanic crisis at Taal volcano (Philippines). *Earth Planets and Space*, 70(1), 1–23. https://doi.org/10.1186/s40623-018-0925-2

Zlotnicki, J., Vargemezis, G., Johnston, M. J. S., Sasai, Y., Reniva, P., & Alanis, P. (2017). Very-low-frequency resistivity, self-potential and ground temperature surveys on Taal volcano (philippines): Implications for future activity. *Journal of Volcanology and Geothermal Research*, 340, 180–197. https://doi.org/10.1016/j.jvolgeores.2017.04.020

## References From the Supporting Information

Chen, C. W., & Zebker, H. A. (2002). Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models. *IEEE Transactions on Geoscience and Remote Sensing*, 40(8), 1709–1719. https://doi.org/10.1109/tgrs.2002.802453

Fattahi, H., & Amelung, F. (2013). DEM error correction in InSAR time series. *IEEE Transactions on Geoscience and Remote Sensing*, 51(7), 4249–4259. https://doi.org/10.1109/tgrs.2012.2227761

Fielding, E. J., Liu, Z., Stephenson, O. L., Zhong, M., Liang, C., Moore, A., & Simons, M. (2020). Surface deformation related to the 2019 M w 7.1 and 6.4 ridgecrest earthquakes in California from GPS, SAR Interferometry, and SAR Pixel Offsets. *Seismological Research Letters*, 91(4), 2035–2046.

Fukuda, J., & Johnson, K. M. (2010). Mixed linear—Non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters. *Geophysical Journal International*, 181(3), 1441–1458.

Jolivet, R., & Agram, P. (2012). Python-based atmospheric phase screen mitigation using atmospheric re-analysis, 2012. Retrieved from http://pyaps.googlecode.com

Rosen, P. A., Hensley, S., Peltzer, G., & Simons, M. (2004). Updated repeat orbit interferometry package released. *Eos, Transactions American Geophysical Union*, 85(5), 47. https://doi.org/10.1029/2004eo050004
<figure><img src="image_10.png" /></figure>

**The 2020 eruption and the large lateral dike emplacement at Taal volcano, Philippines: Insights from radar satellite data**

M. G. Bato¹, P. Lundgren¹, V. Pinel², R. Solidum³, A. Daag³, M. Cahulogan³

¹ Jet Propulsion Laboratory, California Institute of Technology, California, USA

² Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, Grenoble, France

³ Philippine Institute of Volcanology and Seismology (PHIVOLCS), Department of Science and Technology (DOST), Quezon City, Philippines

**Contents of this file**

Text S1

Figures S1 to S17

Table S1

**Introduction**

This supplementary file contains text, figures and a table in support of the main document.

Text S1: Text describing a detailed version of our method section in the main text.

Figures S1-S2: Supporting figures for the pre-eruptive modeling.

Figure S3: Pre- to post-January 2020 eruption volume estimates.

Figure S4-S7: Supporting figures for the co-eruptive modeling.

Figures S8-S12: Supporting figures for the post-main eruptive phase time-series analyses.

Figures S13-S17: Supporting figures for the post-main eruptive phase modeling.

Table S1: Prior conditions used for the inversions as well as the derived MAP solutions from the posterior PDFs.
Text S1. Detailed data processing and modeling

**Pre-eruptive dataset.** For the ALOS-2 SAR data, Earth curvature and topographic effects were corrected based on precise orbits and the SRTM 30-m digital elevation model (DEM). For the Sentinel 1A/B SAR data ARIA uses SRTM 90 m DEM for topographic corrections and relies on precise-orbit ephimerides for orbital correction. All the interferograms were unwrapped using the statistical-cost network-flow algorithm, SNAPHU (Chen & Zebker, 2002) software implemented in ISCE and ARIA. A total of 567 descending and 1004 ascending interferograms were produced by ARIA covering Taal from late 2014 up to January 2020 (i.e. S1 A142: 27 October 2014-11 January 2020, S1 D032: 19 October 2014 – 09 January 2020). For ALOS-2, we produced 12 interferograms for the descending track from 02 December 2018 to 12 January 2020 and nine interferograms for the ascending track from 18 September 2018 to 10 December 2019.

For the Sentinel-1 datasets, we removed interferograms that have an average coherence < 0.6 within the Volcano Island (VI). We also masked pixels that have coherence values that are lower than 0.3-0.35 for both tracks. We applied tropospheric delay corrections using ERA-5 weather models (Jolivet & Agram, 2012) and DEM error corrections (Fattahi & Amelung, 2013) for Sentinel-1, which are all implemented in the MintPy software (Yunjun et al., 2019). For the ALOS-2 data, we did not apply any corrections and masking. The effect of weather models in the pre-eruptive dataset have very minimal impact on the time-series as the amplitude of the surface displacement is larger than the weather-model (and DEM-error) corrections (see Figure S19). However, the corrections were important during the near-real time assessment and monitoring of the volcanic state after the main-eruptive phase.

The NSBAS (Doin et al., 2011) interferograms were 1) produced using the Repeat Orbit Interferometry Package (ROI_PAC) (Rosen et al., 2004) modified to allow TOPSAR data ingestion (Grandin, 2015), 2) unwrapped using the ROI_PAC branch-cut unwrapping algorithm, and 3) corrected for tropospheric delays based on ERA-5 weather models.

**Co-eruptive datasets.** To produce the InSAR phase maps, SAR data from the descending and ascending orbits were processed using ISCE. During the peak of the crisis, the orbital errors were removed using the restituted orbits provided by the Sentinel-1 Quality Control Subsytem to immediately produce the displacement maps and deliver them at minimal latency, ~3 hours, to PHIVOLCS. However, for this manuscript, we reprocessed all the co-eruptive interferograms and use the precise orbits to correct for orbital errors. We unwrapped the phase using SNAPHU and masked the regions where fringe aliasing are observed. Fringe aliasing arises because of the dense number of fringes that often result from large ground displacements.

To generate the pixel offsets map, we applied pixel offset tracking analysis to the Sentinel-1 SAR images on both tracks implemented within ISCE. The pixel offset analysis cross-correlates the co-registered SAR amplitude images to measure the distortion along the radar line-of-sight (slant range) and along-track (azimuth) directions (Fielding et al., 2020). Because of the lower sampling resolution in the azimuth direction (14 m) compared to the slant range (2.3 m) and the near-parallel orientation of the dike to the azimuth direction,
we only used the range offsets for our analysis. Pixel offset maps complement the InSAR
results particularly in regions where fringe aliasing is observed.

**Post-eruptive datasets.** We followed similar processing steps as those of the co-eruptive method using ISCE to produce the interferograms except that we did not apply a mask. We generated 87 and 69 interferograms from the descending and ascending tracks, respectively.

**Pre-eruptive modeling.** We used the Caltech-JPL-developed AlTar v2.0 Bayesian inversion software. AlTar is based on the Cascading Adaptive Transitional Metropolis in Parallel method that combines the Metropolis algorithm with elements of simulated annealing and genetic algorithms, allowing the parallel sampling of high-dimensional problems (Minson et al., 2013). We ran 1000 chains x 1000 steps ($10^6$ samples) starting from uniform distributions (Table S1). Each of the parameter estimates take the form of a *posteriori* probability density functions (PDFs) accounting for the data likelihood and inaccuracies in the data and our prior knowledge about the parameters.

**Co-eruptive modeling.** For the (1) deflating reservoir and a dike with a uniform opening: we performed Bayesian inversion using the Markov chain Monte Carlo (MCMC) method following Fukuda and Johnson (2010) to find the optimal set of parameters that fit our data (i.e. 10 parameters for the CDM, 8 for the rectangular dislocation, 2 for the InSAR shifts). We started from uniform distributions for each of the model parameters (Table S1) and kept $8 \times 10^5$ models from the inversion (i.e. we ran $10^6$ iterations and burned the first 200,000 models). We derived the optimal values by calculating the MAP solutions. For the (2) deflating reservoir and a dike with distributed opening: we solved for a distributed dike model using a non-negative least squares method (Lundgren et al., 2013). We extended the fault width to 20 km and then discretised this into smaller triangular dislocations before applying a Laplacian smoothing operator to regularize the inversion. We tried several smoothing factors to find the balance between the spatial roughness of the tensile opening and the reduction in the root mean square error.
Pre-eruptive Source Modeling

<figure><img src="image_11.png" /><figcaption>**Figure S1.** Histograms of each model parameter after the pre-eruptive source modeling using CDM. X- and Y-loc are the x and y locations of the source in meters relative to 14°N and 121°E. The semi-axes lengths are given by [a, b, c], the rotations are the [ωX, ωY, ωZ], the InSAR offsets for each dataset are [S1, S2, S3, S4] and the opening is *u*. The values in red are the maximum *a* *posteriori* values (MAP) of each model parameters.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure S2.** Pre-eruptive source inversion results. Data (left), model (center), and residual (right) for the quadtreeed ALOS-2 descending track 027, ALOS-2 ascending track 138, Sentinel-1 descending track 032, Sentinel-1 ascending track 142 displacements, arranged from top to bottom. The white star is the location of the pre-eruptive magmatic source. The axes are in meters relative to 14°N and 121°E.</figcaption></figure>
## Volume Estimates

<figure><img src="image_13.png" /></figure>

**Figure S3.** The estimated pre- to post-eruptive volume changes of the CDM (left) and the dike (right) derived from the posterior PDFs of the model parameters. The volume of the CDM (magma reservoir) is calculated using: $ΔV = 4 u (a b + b c + a c)$. Whereas for the dike volume, we use: $ΔV = l w u_3$. We sampled 10000 models from each parameter's posterior distributions to obtain a PDF of the volume change and quantify their respective uncertainties. For the post-eruptive case, we sampled the volume based on the Case-2 modeling approach (i.e. the CDM geometries and locations are fixed given the co-eruptive MAP results) using only the last 50000 MCMC models.
Co-eruptive source modeling

<figure><img src="image_14.png" /></figure>

**Figure S4.** Histograms of each CDM model parameter after the co-eruptive source modeling (i.e. deflating magma reservoir and a dike intrusion). The red lines mark the MAP values of each model parameters whereas the green lines are the mean values. The easting and northing values are given in kilometers relative to 14°N and 121°E. Notice that the negative opening, *u*, falls at the lower bound of our prior. This results from the strong trade-offs between the negative opening and the geometry of the CDM (i.e. semi-axes lengths and the rotations) due to the lack of data constraints within the Volcano Island and in the lake. Hence, in the text, we only reported the volume change of the reservoir.
<figure><img src="image_15.png" /></figure>

**Figure S5.** Histograms of each dike model parameter after the co-eruptive source modeling (i.e. deflating magma reservoir and a dike intrusion). The red lines mark the MAP values of each model parameters whereas the green lines are the mean values. The easting and northing values are given in kilometers relative to 14°N and 121°E.
<figure><img src="image_16.png" /><figcaption>**Figure S6.** Histograms of the shifts in the InSAR measurements. (Left) Descending track 032 co-eruptive data. (Right) Ascending track 142 co-eruptive data. The red lines mark the MAP values of each model parameters whereas the green lines are the mean values.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure S7.** Co-eruptive source inversion results. Data (left), model (center), and residual (right) for the quadtreeed Sentinel-1 descending track 032, Sentinel-1 ascending track 142 InSAR-derived displacements, and Sentinel-1 descending track 032 range-offset, arranged from top to bottom. The axes are in kilometers relative to 14°N and 121°E.</figcaption></figure>
Post-eruptive deformation time-series

<figure><img src="image_18.png" /></figure>

**Figure S8.** Post-main eruptive phase deformation from descending Sentinel-1 time-series between 15 January and 25 June 2020. (Left) Cumulative deformation map in LOS. Red means that the ground moved toward the satellite, blue corresponds to movement away from the satellite. VI: Volcano Island, PRV: Pansipit River Valley, Ref. Pt.: Reference point. (Right) Deformation time-series at selected points within Taal. Orange dots correspond to the raw time-series, blue dots were corrected for atmospheric noise using ERA-5 weather models, and DEM errors.
<figure><img src="image_19.png" /><figcaption>**Figure S9.** The deformation time-series from the Sentinel-1 descending track after applying atmospheric and DEM error corrections, showing the spatial and temporal evolution of the ground between 15 January and 25 June 2020. Black dot is the reference point. All maps are referenced to the first date (15 January).</figcaption></figure>
<figure><img src="image_20.png" /></figure>

**Figure S10.** Post-main eruptive deformation from ascending Sentinel-1 time-series between 17 January and 27 June 2020. (Left) Cumulative deformation map in LOS. Red means that the ground moved toward the satellite, blue corresponds to movement away from the satellite. VI: Volcano Island, PRV: Pansipit River Valley, Ref. Pt.: Reference point. (Right) Deformation time-series at selected points within Taal. Orange dots correspond to the raw time-series, blue dots were corrected for atmospheric noise using ERA-5 weather models, and DEM errors.
<figure><img src="image_21.png" /><figcaption>**Figure S11.** The deformation time-series from the Sentinel-1 ascending track after applying atmospheric and DEM error corrections, showing the spatial and temporal evolution of the ground between 17 January and 27 June 2020. Black dot is the reference point. All maps are referenced to the first date (17 January).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure S12.** The ERA-5 weather models used in correcting the Sentinel-1 ascending time-series. Notice that on 16 April and 03 June 2020, large atmospheric perturbations (i.e. possibly due to heavy rain/humidity) are observed in this region which may have caused the off-trend points in Figure S10.</figcaption></figure>
## Post-eruptive modeling

Case 1: Fixed CDM and dike geometries and locations using co-eruptive MAP solutions

<figure><img src="image_23.png" /></figure>

**Figure S13.** Post-eruptive model for Case 1. (A) Map view and (B) 3D view of the post-eruptive model covering 15 January-04 February 2020. For this model (i.e. post-eruptive model: Case 1), we fixed the geometries and locations of the CDM and the dike using the co-eruptive MAP results (see Table S1). (C) Solution root-mean-square misfit vs. roughness as a function of the smoothing value used for the dike. The value in red text is the smoothing factor that we used for the post-eruptive model.
<figure><img src="image_24.png" /></figure>

**Figure S14.** Post-eruptive source inversion results for Case 1 (Figure S13). The geometries and locations of the CDM and the dike were fixed using the co-eruptive MAP values presented in Table S1. Afterwards, we inferred for the negative opening of the CDM and the distributed opening of the dike using non-negative least squares approach. Data (left), model (center), and residual (right) for the quadtreeed Sentinel-1 descending track 032, Sentinel-1 ascending track 142 InSAR-derived displacements arranged from top to bottom. The axes are in kilometers relative to 14°N and 121°E.
Case 2: Fixed CDM geometries and locations using co-eruptive MAP solutions. CDM opening and dike parameters are estimated using MCMC

<figure><img src="image_25.png" /></figure>

**Figure S15.** Histograms of each CDM model parameter after the post-main eruptive phase source modeling using MCMC (i.e. deflating magma reservoir and a dike intrusion). Notice that we only inferred for the negative opening of the CDM. The red lines mark the MAP values of each model parameters whereas the green lines are the mean values. The easting and northing values are given in kilometers relative to 14°N and 121°E.
<figure><img src="image_26.png" /></figure>

**Figure S16.** Histograms of each dike model parameter after the post-main eruptive phase source modeling using MCMC (i.e. deflating magma reservoir and a dike intrusion). The red lines mark the MAP values of each model parameters whereas the green lines are the mean values. The easting and northing values are given in kilometers relative to 14°N and 121°E.
<figure><img src="image_27.png" /></figure>

**Figure S17.** Post-eruptive model for Case 2. (A) Map view and (B) 3D view of the model covering 15 January-04 February 2020. (C) Solution root-mean-square misfit vs. fault slip roughness as a function of the smoothing value used for the dike. The value in red text is the smoothing factor that we used for the post-eruptive model.
<figure><img src="image_28.png" /></figure>

**Figure S18.** Post-eruptive source inversion result for Case 2. Data (left), model (center), and residual (right) for the quadtreeed Sentinel-1 descending track 032, Sentinel-1 ascending track 142 InSAR-derived displacements arranged from top to bottom. The axes are in kilometers relative to 14°N and 121°E.
<figure><img src="image_29.png" /><figcaption>**Figure S19.** Displacement time-series at 14.0188 °N, 121.0012 °E of ERA5+DEM correction values (blue) and the corrected time-series (orange, as shown in Figure 2C).</figcaption></figure>
<table><caption>**Table S1.** The prior conditions used during the inversions. In the co- and post-eruptive cases, we first inferred for the geometries and locations of the CDM and the dike using MCMC so that we can fix them before computing for the distributed slip. Here, we present the resulting maximum a posteriori (MAP) values after the MCMC inversions. Note that for the post-eruptive modeling, we either used the co-eruptive MAP values or we performed an independent MCMC inversion for the geometries and locations of the dike.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Pre-eruptive</th><th colspan="2">Co-eruptive</th><th colspan="2">Post-eruptive</th></tr><tr><th>Prior distribution</th><th>Prior distribution</th><th>MAP from MCMC</th><th>Prior distribution</th><th>MAP from MCMC</th><th></th></tr></thead><tbody><tr><td><strong>CDM</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>X-location, X-loc (km)</td><td>U[-10, 10]</td><td>U[-5, 5]</td><td>0.215</td><td>U[-5, 5]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Y-location, Y-loc (km)</td><td>U[-10, 10]</td><td>U[-5, 5]</td><td>1.833</td><td>U[-5, 5]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Depth (km)</td><td>U[1, 10]</td><td>U[0.1, 20]</td><td>5.703</td><td>U[0.1, 20]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Semi-axes lengths:</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a (km)</td><td>U[1, 10]</td><td>U[0.1, 10]</td><td>0.728</td><td>U[0.1, 10]</td><td>Coerupt MAP</td><td></td></tr><tr><td>b (km)</td><td>U[1, 10]</td><td>U[0.1, 10]</td><td>2.312</td><td>U[0.1, 10]</td><td>Coerupt MAP</td><td></td></tr><tr><td>c (km)</td><td>U[1, 10]</td><td>U[0.1, 10]</td><td>0.318</td><td>U[0.1, 10]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Rotations:</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>X-rotation, ω<sub>X</sub> (°)</td><td>U[-30, 30]</td><td>U[-89, 89]</td><td>1.15</td><td>U[-89, 89]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Y-rotation, ω<sub>Y</sub> (°)</td><td>U[-89, 0]</td><td>U[-89, 89]</td><td>-15.71</td><td>U[-89, 89]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Z-rotation, ω<sub>Z</sub> (°)</td><td>U[-89, 0]</td><td>U[-180, 180]</td><td>23.96</td><td>U[-180, 180]</td><td>Coerupt MAP</td><td></td></tr><tr><td>Opening, u (m)</td><td>U[0, 10]</td><td>U[-50, 0]</td><td>-50</td><td>U[-50, 0]</td><td>-0.725</td><td></td></tr><tr><td><strong>Dike</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Length, l (km)</td><td>--</td><td>U[0.1, 40]</td><td>20.966</td><td>U[0.1, 40]</td><td>18.152</td><td></td></tr><tr><td>Width, w (km)</td><td>--</td><td>U[0.1, 20]</td><td>8.152</td><td>U[0.1, 20]</td><td>11.738</td><td></td></tr><tr><td>Depth (km)</td><td>--</td><td>U[0, 15]</td><td>0.556</td><td>U[0, 15]</td><td>0.712</td><td></td></tr><tr><td>Dip (°)</td><td>--</td><td>U[-260, 90]</td><td>-76.8</td><td>U[-260, 90]</td><td>-92.1</td><td></td></tr><tr><td>Strike (°)</td><td>--</td><td>U[-90, 90]</td><td>33.1</td><td>U[-90, 90]</td><td>47.5</td><td></td></tr><tr><td>X-location, Easting (km)</td><td>--</td><td>U[-20, 40]</td><td>-3.512</td><td>U[-20, 40]</td><td>-6.376</td><td></td></tr><tr><td>Y-location, Northing (km)</td><td>--</td><td>U[-20, 10]</td><td>-3.899</td><td>U[-20, 10]</td><td>-5.409</td><td></td></tr><tr><td>Slip, u<sub>3</sub> (m)</td><td>--</td><td>U[-5, 15]</td><td>3.2</td><td>U[-5, 15]</td><td>0.264</td><td></td></tr><tr><td>InSAR offsets (m)</td><td>U[-1, 1]</td><td>U[-1, 1]</td><td>-0.0024 (des);<br>0.0648 (asc)</td><td>U[-1, 1]</td><td>-0.0187 (des);<br>0.0081 (asc)</td><td></td></tr></tbody></table>
References:

Chen, C. W., & Zebker, H. A. (2002). Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models. IEEE Transactions on Geoscience and Remote Sensing, 40(8), 1709-1719.

Doin, M.-P., Guillaso, S., Jolivet, R., Lasserre, C., Lodge, F., Ducret, G., & Grandin, R. (2011). Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using Envisat data. In Proceedings of the fringe symposium (pp. 3434-3437).

Fattahi, H., & Amelung, F. (2013). DEM error correction in InSAR time series. IEEE Transactions on Geoscience and Remote Sensing, 51 (7), 4249-4259.

Fielding, E. J., Liu, Z., Stephenson, O. L., Zhong, M., Liang, C., Moore, A., ... Simons, M. (2020). Surface Deformation Related to the 2019 M w 7.1 and 6.4 Ridgecrest Earthquakes in California from GPS, SAR Interferometry, and SAR Pixel Offsets. Seismological Research Letters.

Fukuda, J., & Johnson, K. M. (2010). Mixed linear —non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters. Geophysical Journal International, 181 (3), 1441-1458.

Grandin, R. (2015). Interferometric processing of SLC Sentinel-1 TOPS data.

Jolivet, R., & Agram, P. (2012). Python-based atmospheric phase screen mitigation using atmospheric re-analysis, 2012. URL: http://pyaps.googlecode.com

Lundgren, P., Poland, M., Miklius, A., Orr, T., Yun, S.-H., Fielding, E., ... others (2013). Evolution of dike opening during the March 2011 Kamoamoa fissure eruption, Kilauea Volcano, Hawaii. Journal of Geophysical Research: Solid Earth, 118 (3), 897-914.

Minson, S., Simons, M., & Beck, J. (2013). Bayesian inversion for finite fault earthquake source models |Theory and algorithm. Geophysical Journal International, 194 (3), 1701-1726.

Rosen, P. A., Hensley, S., Peltzer, G., & Simons, M. (2004). Updated repeat orbit interferometry package released. Eos, Transactions American Geophysical Union, 85 (5), 47.

Yunjun, Z., Fattahi, H., & Amelung, F. (2019). Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction. Computers & Geosciences, 133, 104331.