Geophysical Responses to an Environmentally‐Boosted Volcanic Unrest

The spatiotemporal relationship between geophysical, environmental, and geochemical responses during volcanic unrest is essentially unknown, making their joint use and interpretation for eruption forecasting challenging. Here, Empirical Orthogonal Functions analysis applied to GPS data allows the separation of the dominant deep‐sourced inflation from environmentally controlled signals associated with extension at Campi Flegrei caldera. This separation bridges the gap between deformation, seismic and geochemical responses, clarifying the processes underlying the ongoing volcanic unrest. Persistent meteoric forcing during the 2017–2018 hydrological year changed the decadal trend of seismic energy and secondary deformation components, pairing their spatial patterns. The result was a block in the carbon dioxide released in 2018 at Solfatara, the primary stress‐release valve at the caldera. The subsequent overpressure weakened the fractured eastern caldera, opening pathways for deep, hot materials to reach the surface. Our results give insight into how environmental forcing can favor volcanic unrest in pressurized calderas.


Introduction
The method of Empirical Orthogonal Functions (MEOF), a powerful tool for weather forecasting (Lorenz, 1956) and climate and oceanic modeling (Hannachi et al., 2007;Navarra & Simoncini, 2010;Smith et al., 1996), is equivalent to a Principal Component Analysis applied to data that have a combination of spatial and temporal trends (Zerbini et al., 2013).MEOF decomposes a space-time field into spatial patterns and associated time indices, where the components reflect different but interconnected processes (Navarra & Simoncini, 2010).When applied to GPS data at the continental scale over a decade, the main EOF components reveal environmental forcing associated with a hot, dry season (Zerbini et al., 2013).Its application to InSAR data at Campi Flegrei caldera (Southern Italy, Figures 1a and 1b) shows that the first EOF component captures most of the deformation signal and can be modeled as the response to magma or magmatic fluid inflation, producing unrest at least since 2015 (Amoruso & Crescentini, 2022).In a volcano, transient magmatic, geochemical, and environmental processes are tightly entangled; it is thus essential to analyze components other than the first, as they might have a higher sensitivity to such transients (Monahan et al., 2009).
Forcing similar to that observed at the continental scale (Zerbini et al., 2013) was most likely active during the 2017-2018 wet hydrological year (Table S1 in Supporting Information S1) when anomalous continuous precipitations hit the caldera during the strongest drought of the decade (2015-2020-Figure 1c).Deeper magmatic reservoirs are primary contributors to seismicity at Campi Flegrei; however, seismic energy (Tramelli et al., 2021; Figure 1b and Figure S1 in Supporting Information S1) steepens after July 2014 (from gray to black seismic energy trend) and more clearly in September 2017 (from black to red, Figure S1 in Supporting Information S1).Short-range (from hours to days) seismic rates and intensities are tied to external modulators like tides and precipitations (Hainzl et al., 2006;Matthews et al., 2009;Neuberg, 2000) in such pressurized hydro-geothermal systems (Petrosino & De Siena, 2021;Scafetta & Mazzarella, 2021).Carbon dioxide (CO 2 ) emissions at Solfatara and its faults (Bruno et al., 2007;Di Martino et al., 2022, S-Figures 1b, 2a, and 3a) are primary manifestations of the Campi Flegrei unrest.The second portion of the 2017-2018 hydrological year (February-September) also corresponds to the longest and largest CO 2 flux sink ever recorded at Solfatara (Figures 1b and 1c and Figure S1 in Supporting Information S1).
The relationship between deformation and rainfall is more ambiguous (Petrosino et al., 2021).Raw GPS signals show variations that are difficult to link to coeval geochemical and seismic signals across the whole caldera and are primarily associated with deep-sourced carbon or magma throughout the 2005-2019 unrest (Amoruso & Crescentini, 2022;Amoruso et al., 2015;Moretti et al., 2018;Vanorio & Kanitpanyacharoen, 2015).The injection of meteoric fluids likely leads to thermo-poroelastic processes (Moretti et al., 2018;Nespoli et al., 2021;Scafetta & Mazzarella, 2021), but separating the deformation signals produced by deep and shallow processes and linking them to seismicity remains challenging.However, the anomalous continuous precipitations in the middle of the 2015-2020 drought (Figure 1c) represent a long-range environmental forcing that the MEOF can detect (Zerbini  (Isaia et al., 2015;Vilardo et al., 2007) is superimposed on the yellow and orange isosurfaces representing high-V p /V s (1983-1984-Vanorio et al., 2005) and low-velocity (2011-2013-De Siena et al., 2018) anomalies, respectively.The bottom of the high-attenuation isosurface (gray, 1983-1984) coincides with a hot inflating and degassing reservoir (Buono et al., 2022;Castaldo et al., 2021).Green symbols (+, Pepe et al., 2019): seismic front of the thermo-poroelastic reservoir in 2012-2014.Colored circles: earthquakes above magnitude 0.75 color-coded by nucleation time.P = Pozzuoli; S = Solfatara; Pi = Pisciarelli; M = Monte Nuovo; A = Astroni; G = Gauro.Arrows at the surface: GPS displacement directions and total uplift on 31 December 2019 (De Martino et al., 2021).(b) From top to bottom, the temporal relation between soil CO 2 flux at Solfatara (Chiodini et al., 2021), EC1, EC2, and Z-score (seismic energy, Tramelli et al., 2021) from 22 July 2014 to December 2019.Changes from gray to black to red in the time series correspond to three different seismic energy trends associated with varying seismic swarm durations (see Figure S1 in Supporting Information S1).Yellow bars and dates highlight the 2018 sink (1 March-29 September).Thin black and orange lines mark 01/03 and 21/06 every year, respectively.(c) From top to bottom, the temporal relation between average daily precipitations over a month, and EC2 and EC3 calculated from weekly GPS data extending to December 2022 (De Martino et al., 2023).Black vertical lines represent the start of a hydrological year.The blue-colored lines show trends from 1 September 2020, onwards.Dry (D), Very Dry (VD), and Wet (W) years are reported with precipitations and defined following Table S1 in Supporting Information S1.  et al., 2013) once applied to daily GPS data within the caldera (Figure 1a, GPS measurements, De Martino et al., 2021).
This study uses the MEOF to identify deep-sourced and environmental processes underlying the GPS deformation field.It shows that secondary deformation patterns linked to environmental forcing pair with seismic migrations, while their temporal components detect the closure of the primary stress-relief valve at Campi Flegrei: outgassing at Solfatara.Because of the subsequent shallow overpressure, a long-lived front (Pepe et al., 2019) expanded east, widening the fumarolic field at Pisciarelli (Fedele et al., 2020) and breaking the weakened crust in 2019, boosting the ongoing deep-sourced unrest.

Empirical Orthogonal Functions
GPS displacements (North, East, and Up) measured across the caldera are arranged in a data matrix X having as many columns as the number of GPS stations multiplied by the number of the displacement components (three) and as many rows as the number t of measurement times.We limit our analysis to co-recording stations having no extended gaps (Figure S2a in Supporting Information S1 shows stations used from 2013/2014 and 2016).By singular value decomposition, we obtain: (1) Navarra & Simoncini, 2010) give the temporal patterns; the columns of the matrix EOF give the (mutually orthogonal) spatial patterns.We thus associate a temporal (ECn) and a spatial (EOFn) pattern to each component n at each station.The relevance of both patterns depends on the diagonal elements in matrix D or (as more commonly presented) on Λ k = λ 2 k , the covariance matrix (Figure S2b in Supporting Information S1).We estimate a component significance from the eigenvalue errors, given as Λ k sqrt(2/p), where p is the number of statistically independent samples (Navarra & Simoncini, 2010) calculated using the autocorrelation function (Thiiébaux & Zwiers, 1984).Usually, few components capture most of the original information; however, the need for spatial orthogonality, temporal uncorrelation, and nonlocal maximization of variance over the domain in different components implies that each component has no dynamical, kinematic, and statistical meaning independently of the others (Monahan et al., 2009).
A quasi-linear trend induced by the Campi Flegrei uplift dominates the GPS time series; thus, we first compute the autocorrelation function of X to estimate p for EC1; then, we estimate p for the other components combined (i.e., ECN) by computing the autocorrelation for a matrix obtained by subtracting EC1 from X.We noted an unrealistic correlation between specific temporal variations of EC1 and EC2 at all stations of the caldera (Figure S4 in Supporting Information S1).This effect is probably due to the slight high-frequency movement (noise) of the reference GPS stations.We thus analyzed displacement-removed data: we removed the mean of each displacement component over n r stations outside the caldera from each row in X (i.e., at each time).The displacements of the external stations are arranged in a data matrix R having n t rows and 3n r columns.For each X i,j , we compute k = j%3 (% is the modulo operator, and k = 3 if j%3 = 0).The elements of the displacement-removed data matrix Y are: (2)

Earthquake Locations
The Campi Flegrei seismicity is monitored by the permanent and mobile seismic networks of Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Napoli-Osservatorio Vesuviano (INGV-OV).The data set consists of volcano-tectonic earthquakes detected in 2005-2019 (color-coded circles, Figure 1a).P-and S-phases recorded by the networks were picked, integrated, and inverted with a probabilistic non-linear 3D grid search.We used the NLLOC software (Lomax et al., 2014) to search for the best solution within the 3D velocity structure of Battaglia et al. (2008).Locations with RMS higher than 0.1 were removed from the data set.

Deformation Modes and Joint Spatial Analysis
We  S2a, arrows in Supporting Information S1), analogue to the one obtained using InSAR data (Amoruso & Crescentini, 2022).The dependence of GPS data on inflation means that few samples are independent when including EC1, which is always statistically significant (Figure S2b in Supporting Information S1).The autocorrelation computed on ECN shows that only EC2 is statistically significant (insets in Figure S2b in Supporting Information S1).The displacement-removed processing drastically reduces noise on EC1 and its correlations with EC2 (Figure S5 in Supporting Information S1).
From now on, we will always refer to displacement-removed components.EC2 is our primary proxy for the time variation of components different from EC1 due to its statistical significance and contribution to the variance (Figure S2b in Supporting Information S1).It shows similar oscillatory behaviors in all periods considered, independently from differences in station numbers (Figure 1b, Figure S2c in Supporting Information S1).EOF1 shows inflation with a maximum at the center of the 1983-1984 reservoir in all periods (Figure S2a in Supporting Information S1).It is almost identical to the GPS pattern (Figure 1a).From 2015, EC2 shows an annual cycle (Figure 1b) comprising a decrease from the start of March (black vertical lines) to 21 June (orange lines), a plateau through the summer, and an increase from September to the end of February (black and red curved arrows in Figure 1b).The only years with no plateau are 2018 and 2019, with a decrease lasting ∼1 and ∼2 months longer (Figure 1b) followed by a sudden increase during the sink.We also show EC2 calculated from 2016 with seismic energy around the period of CO 2 sink in 2018 (Figure 2a) and during the subsequent seismic unrest (Figure 3a).yellow = youngest).A seismic front appears when seismicity consistently distributes across a volume without an apparent temporal trend.

Comparison Between Temporal Trends
Scafetta and Mazzarella (2021) show that a 7-11 September 2017 extreme rain event (>200 mm of meteoric water over 4 days, twice the total amount of an average September) follows a ∼2-year gap in the correlation between seismic swarms and rainfall, corresponding to a 2-year drought (Figure 1c).At the beginning of 2019, the CO 2 flux recovered the exponential trend followed before the sink.EC1 and EC2 show no apparent relation to these rainfall events, with our processing removing the minor, transient increases at both ends of the sink (Figure S5 in Supporting Information S1, EC1 and EC2).The 2018 cycle coincides with the beginning of the CO 2 sink (1 March 2018, Figures 1c and 2a).The 2019 seasonal EC2 decrease comprises the maximum CO 2 flux measured at Solfatara until 2019 (Figure 1c); if calculated from 2016, the 2019 decrease ended in August (Figure 3a), ∼2 months later than in 2016 and 2017, which suggests a primary change in the behavior of the EC2 time series hidden by its seasonal oscillations.

Deformation and Seismic Responses to a Carbon Sink
Deep deformation sources have been active since at least 2015 at the caldera and are controlling EC1 (Amoruso & Crescentini, 2022).The September 2017 events were the first manifestation of some of the highest (Table S1 in Supporting Information S1) and most continuous yearly precipitations recorded since 2012 at the caldera (Figure 1c).They happened in the middle of a 5-year drought, which, at the continental scale, controls EOF components (Zerbini et al., 2013).EC2 thus appears to be a significant marker of year-long environmental and seasonal processes alternative to deep-sourced deformation.However, it must be considered with additional components to infer their dynamical, kinematic, and statistical meaning.Hence, we compare their combined spatial (EOFN-arrows) and temporal (the ECN color scale shows the average displacement over the period) patterns with coeval seismic migrations and fronts in Figures 2b-2d and 3b-3d (see ECN and EOFN across the entire caldera in Figure S3 in Supporting Information S1).
In the 5 months before the sink, earthquakes consistently migrated from NE (Pisciarelli) to SW across the SW-to-NE-striking fault crossing Solfatara (color-coded arrow, Figure 2b).Seismic activity on this fault has been a primary stress-release valve at the caldera since the unrest started, in 2005 (Bruno et al., 2007;Chiodini et al., 2021;Isaia et al., 2015;Vilardo et al., 2007).The fault remains the only seismically active, degassing connection between the expanding central reservoirs and a critically stressed front, contouring the eastern reservoir boundary between 2012 and 2014 (Pepe et al., 2019-Figure 2d, green + symbols).Right before the sink, these migrations were associated with compression of the 1983-1984 reservoir (ECN sinks into the reservoir, Figure 2b), while the stress from deeper sources increased at a higher rate (EC1, Figure 1c).The EC2 pattern points to a compression of the 1983-1984 supercritical fluid reservoir, which coupled with the extensional mechanics (Tramelli et al., 2022) responding to inflation and deflation from the deep, hot reservoir (Amoruso & Crescentini, 2022).
The shallow (∼1 km depth) volumes under the Solfatara fault remain aseismic throughout the sink and right after it (Figures S6c and S6d in Supporting Information S1, red S).Solfatara presents a clay cap with few faults, allowing fluid migration (Bruno et al., 2007;Isaia et al., 2015;Siniscalchi et al., 2019;Vilardo et al., 2007).Increased pore pressure from the 2017-2018 precipitations, following the much dryer 2016-2017, can increase permeability across the cap by opening cracks and pores; however, progressive increases in pore pressure dislodge clay minerals and impede the flow from depth (e.g., Coyner, 1977) so that the cap becomes a plug in the carbon flow (Figure 2a), increasing stress as the system dries up in the summer (Moretti et al., 2018).This process is enhanced by thermal expansion and liquefaction (Chiodini et al., 2021;Siniscalchi et al., 2019), which further lower the ability of rocks to fracture, producing earthquakes.While the relation between gas emissions at Solfatara and Pisciarelli was unclear before the sink, their time series were anticorrelated during the sink (Figure S8 in Supporting Information S1).The plug caused earthquakes to plateau at ∼1 km under Solfatara while seismic fluids migrated toward Pisciarelli (Figure 2c) through communicating aquifers (Isaia et al., 2022).During the sink, CO 2 -bearing fluids could only surface near the Pisciarelli vent, where the clay cap had not formed (Siniscalchi et al., 2019), drastically increasing the corresponding emissions.

Geophysical Research Letters
10.1029/2023GL104895 DE SIENA ET AL.
The relationship between EOFN, seismic patterns, and reservoir locations clarifies how the system released the combined stress from the plugged Solfatara and the deep inflating sources at the end of the sink.During the sink, seismic migrations rose from Pisciarelli.They contoured the 1983-1984 reservoir (Figure 2c, splitting colorcoded arrows, maximum earthquake depths at 2.3 km), defining its eastern 2012-2014 front (Nespoli et al., 2021;Pepe et al., 2019).This seismic behavior is typical of stress release from an over-pressurized thermoporo-elastic reservoir, specifically of the opening of its shallow crustal front (Lu, 2018;Nespoli et al., 2021).The 2011-2013 low-velocity reservoir was depressurized in 2018 (Petrosino & De Siena, 2021), storing fluids that migrated toward the east.On the contrary, the 1983-1984 supercritical fluid reservoir remained aseismic and coincided with the area of highest pressurization throughout the unrest (Petrosino & De Siena, 2021).
The result of the opening at the top of the system was the definition of a new seismically active front at the end of 2018, further north and east of the central supercritical fluid reservoirs (Figure 2d, green curve), which released stress and stopped the CO 2 flux sink (Figure 2a).Right after the sink, seismic patterns distributed much further east and north than the pre-existing front (+symbols, Figure S6d in Supporting Information S1), almost reaching the Astroni crater.In the same period, the EOFN started rotating toward the east with the corresponding ECN (color scale) lower than during the sink (see Figures 2b-2d), a signal of the incoming predominance of deepsourced processes.The shift of the seismic front and the start of a clockwise rotation of the ECN are a symptom of the incoming extensional dynamics that will produce the highest-magnitude earthquakes at the caldera since 1984 (Tramelli et al., 2022).

Top-Down Control on Stress Release Contributes to Unrest
The opening of the shallow front was a vital contributor to the 2019 unrest, characterized by excessive outgassing (Figure 3a) and the highest seismic and deformation rates recorded since 1984 (Chiodini et al., 2021).Recovery from the sink corresponded with the re-establishment of seismic migrations from the caldera center to Pisciarelli across the Solfatara fault, coeval to the release of extensional stress from the 1983-1984 reservoir (Figure 3b, color-coded arrow).From this point onwards, EC2 and seismic energy followed opposite linear trends (Figure 3a) until the 3.1 Md event of 6 December 2019.The ECN spatial patterns had fully rotated east (Figure 3b), imaging expansion of the central reservoir associated with co-directed, downward seismic migrations.We infer that ECN and seismic migrations are proxies for extensional dynamics, with their patterns opposite to those at the start of the sink.They show (a) the close link between secondary EC components and the release of seismic energy and (b) how opposite geophysical responses characterize geochemical sink and recovery.
The Solfatara plugging and associated liquefaction and thermo-elastic processes (Scafetta & Mazzarella, 2021;Siniscalchi et al., 2019) had weakened the shallow portion of the caldera east of the supercritical reservoir in 2018, but seismic migrations had only cracked the eastern boundary of the shallow reservoir and opened a shallow front.
In 2019, the result of this shallow caldera staging (Roman & Cashman, 2018) was (a) the rise of materials from the deep degassing reservoir to the eastern tip of the main SW-NE striking fault at Solfatara and (b) the seismic opening of the SN pressurized transfer structure connecting Pozzuoli and Astroni and acting as new seismic front (Figure 3b, green curve- Petrosino & De Siena, 2021).The recovery of the CO 2 flux coincided with consistent seismic migrations east of the front opened at least since 2014 (+symbols, Figure S9 in Supporting Information S1), which will generate the highest-magnitude earthquakes since 1984 at the caldera.ECN patterns in this period pointed to the deep, hot source of unrest right below Pozzuoli, opening a deeper seismic front (Figure 3c) in the eastern caldera and beyond.The mechanisms acting before 2017-2018 could not sustain the corresponding stress, leading to higher seismic rates in 2019 (Tramelli et al., 2021).The EC1 rates are the strongest across our data set.In this period, the new seismic front expanded east, while EOFN points to the volumes associated with the highest-magnitude earthquakes (Figures 3d, Tramelli et al., 2021).
On 8 November 2023, weekly GPS data extending to November 2022 were released (De Martino et al., 2023), allowing us to process data over the same array of stations, even if at a lower time resolution.The MOEF effectively separates the long-range environmental and extensional behavior (EC2, Figure 1c) previously hidden under the oscillatory short-range seasonal forcing (EC3, Figure 1c): these two components are the primary contributors to the EOFN patterns, spatially coupled with seismic migrations and fronts (Figures 2 and 3).The two leading secondary EOF components (Figure S10 in Supporting Information S1) thus depend on environmental forcing, with the EC2 stress release (decrease) starting after the sink at Solfatara, that is, at the end of the 2017-2018 wet hydrological year (Figure 1c).This is when EOFN patterns rotate while seismicity marks a shallow front expansion (Figures 2c and 3a), the turning point for deep materials to rise to the surface (Figures 3b-3d).A stronger decrease characterizes EC2 around the start of the wet 2020-2021 hydrological year (Figure 1c, blue), but the lower temporal resolution of GPS data forbids a joint spatial analysis with seismicity.Our results suggest it is possible to detect changes leading to seismic unrest in the eastern caldera from geochemical (CO 2 sink, Figure 1b) and secondary deformation (Figure 1c, EC2) time series.Analyzing decadal daily GPS recordings could lead to more accurate forecasts and clarify the contribution of seasonal and decadal environmental effects on deformation data.

Conclusions
This study represents the first joint analysis of Empirical Orthogonal Functions applied to GPS data, time-lapse earthquake migrations and CO 2 flux time series.Its results show the geophysical time-resolved response to environmental and deep processes producing unrest at a volcanic caldera-an image of how the volcano responds to being increasingly stressed from depth and surface.The caldera momentarily reduces its carbon outflow, storing rising carbon-bearing fluids within its eastern faults.This storage becomes an additional pressure source that can further crack the shallow fractured eastern caldera around expanding thermo-poroelastic supercritical fluid reservoirs, opening a pathway for deeper materials to rise to the surface.
The relationships between the 1983-1984 supercritical fluid reservoir, the ECN and EOFN secondary deformation patterns, and seismic migrations show that extensional processes associated with fluid migrations are a vital release valve for increases in thermo-poroelastic stress primarily caused by deep degassing sources.Increased shallow carbon storage and year-long environmental forcing can perturb this cyclic stress release mechanism.Our study suggests that continuous long-lasting meteoric forcing amid extended periods of drought contributes to geophysical and geochemical responses primarily controlled by deep processes: understanding this forcing appears central to assessing volcanic hazards at this and other volcanic calderas.at the INGV-OV for providing continuous routine seismic peaks and the 3D probabilistic locations of the Campi Flegrei earthquakes.We are grateful to Antonio Navarra, Prospero De Martino and an anonymous reviewer for providing data and comments that allowed us to strengthen our results.Open Access funding enabled and organized by Projekt DEAL.

Figure 1 .
Figure 1.Geomorphology, seismicity, deformation, and geochemistry at Campi Flegrei.(a) The geomorphological and structural map of the caldera(Isaia et al., 2015;Vilardo et al., 2007) is superimposed on the yellow and orange isosurfaces representing high-V p /V s(1983-1984-Vanorio et al., 2005) and low-velocity (2011-2013- De Siena et al., 2018) anomalies, respectively.The bottom of the high-attenuation isosurface(gray, 1983-1984)  coincides with a hot inflating and degassing reservoir(Buono et al., 2022;Castaldo et al., 2021).Green symbols (+,Pepe et al., 2019): seismic front of the thermo-poroelastic reservoir in 2012-2014.Colored circles: earthquakes above magnitude 0.75 color-coded by nucleation time.P = Pozzuoli; S = Solfatara; Pi = Pisciarelli; M = Monte Nuovo; A = Astroni; G = Gauro.Arrows at the surface: GPS displacement directions and total uplift on 31 December 2019 (De Martino et al., 2021).(b) From top to bottom, the temporal relation between soil CO 2 flux at Solfatara (Chiodini et al., 2021), EC1, EC2, and Z-score (seismic energy, Tramelli et al., 2021) from 22 July 2014 to December 2019.Changes from gray to black to red in the time series correspond to three different seismic energy trends associated with varying seismic swarm durations (see Figure S1 in Supporting Information S1).Yellow bars and dates highlight the 2018 sink (1 March-29 September).Thin black and orange lines mark 01/03 and 21/06 every year, respectively.(c)From top to bottom, the temporal relation between average daily precipitations over a month, and EC2 and EC3 calculated from weekly GPS data extending to December 2022(De Martino et al., 2023).Black vertical lines represent the start of a hydrological year.The blue-colored lines show trends from 1 September 2020, onwards.Dry (D), Very Dry (VD), and Wet (W) years are reported with precipitations and defined following TableS1in Supporting Information S1.

Figure 2 .
Figure 2. Seismic and secondary deformation patterns before and during the sink.(a) CO 2 soil flux, EC2, and seismic energy time series from the 2017 rainfall events (Scafetta & Mazzarella, 2021) to the end of the sink.Line colors are identical to those in Figure 1b.Green color scale: ECN average deformation in each period.(b-d) Comparison between the EOFN and ECN patterns (arrows at the surface and their color), interpreted seismic migrations (arrows color-coded like earthquakes in Figures S6b-S6d in Supporting Information S1) and fronts.Solfatara (S) is labeled in red during sink periods and otherwise black.
apply the MEOF to deformation signals recorded until December 2019: (a) from January 2013, (b) from July 2014, and (c) from July 2016 (see Text S1 in Supporting Information S1).In all cases and independently of processing: (a) EC1 explains >98.9% of the variance (Figure S2b in Supporting Information S1); (b) EOF1 patterns correspond to the inflation of deformation sources in the center of the caldera (Figure

FiguresFigure 3 .
Figures 2b-2d and 3b-3d show jointly: 1.The summed spatial pattern of the displacement of all components except the first (EOFN, arrows, also named Extensional/Reservoir component) for stations inside the caldera rim.The associated green and white color scale shows the average ECN displacement over the same period.Figures S3b-S3g in Supporting Information S1 shows ECN and EOFN at all stations used in the same periods.2. Earthquake migrations (color-coded arrow) and fronts (green curves) defined from earthquake locations (circles, Figures S6 and S7 in Supporting Information S1) spanning the same period.In Figures 2 and 3, locations and interpreted migrations are labeled using the same color scale (cyan = oldest;