# Modeling magma recharge dynamics during the 2016 Nevados de Chillan eruption: An interacting two-chamber system evidenced by petrology and geodesy

Camila Novoa Lizama, D. Remy, J.C. Baez, A. Oyarzun, Sylvain Bonvalot, A.

Hooper

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

## ► To cite this version:

Camila Novoa Lizama, D. Remy, J.C. Baez, A. Oyarzun, Sylvain Bonvalot, et al.. Modeling magma recharge dynamics during the 2016 Nevados de Chillan eruption: An interacting two-chamber system evidenced by petrology and geodesy. Journal of Volcanology and Geothermal Research, 2025, 458, pp.108253. 10.1016/j.jvolgeores.2024.108253. hal-05192546

HAL Id: hal-05192546

https://hal.science/hal-05192546v1

Submitted on 30 Jul 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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Modeling magma recharge dynamics during the 2016 Nevados de Chillan eruption: An interacting two-chamber system evidenced by petrology and geodesy

Camila Novoa Lizama<sup>a,*</sup>, D. Remy<sup>b</sup>, J.C. Baez<sup>c</sup>, A. Oyarzun<sup>d</sup>, S. Bonvalot<sup>b</sup>, A. Hooper<sup>a</sup>

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

## ARTICLE INFO

**Keywords:**

InSAR and GNSS

Explosive and effusive eruption phases

Southern Volcanic Zone

Double reservoir modeling

Uplift episode

Phreatomagmatic eruption

## ABSTRACT

The Nevados de Chillan is a large composite stratovolcanic complex located in the Southern Andean Volcanic Zone, at the south of Chile. Its last eruption spanned six years from January 8, 2016 until January 2023. After three years without any deformation, exhibiting phreatic and phreatomagmatic activity, an uplift episode was reported in June 2019, marking the start of its magmatic phase. How this geodetic activity evolved in the following three years is still not clear, neither how it correlates with its superficial activity. Here we analyze InSAR time series, as well as daily GNSS time series from 2015 to early 2022, highlighting the wide range of ground surface displacements observed during the eruption. In-depth analysis of these displacements together with recent petrological and geochemical results leads us to consider a recharge mechanism involving a double-reservoir model to explain the observed geodetic activity. We develop an analytical model of dynamic magma flow, which we coupled with a boundary element method to account for any geometry of the reservoirs and the topography of the volcano. The model consists of a shallow elongated source located at 5.8 km depth below the volcanic edifice connected by an incompressible magma-filled hydraulic pipe to a deeper sill like source centered at 15 km depth. We propose that the activation of the system started with a small magma intrusion of one month of duration, which was sufficient to overheat the hydrothermal system and re-mobilize magma in the shallow chamber, thus explaining the non-deformation during its phreatic phase and the slight subsidence observed during its phreato-magmatic stage. Then, we propose that a new and larger magma intrusion occurred in June 2019, explaining the uplift episode observed, which continued for the following three years of the eruption, decaying exponentially. Our model indicates that this intrusion was triggered by magma coming from the crust-mantle boundary to the deep reservoir at constant rate of 0.016 km³y⁻¹ from June 2019 to Jan. 2022, with small changes to this rate that would explain the small fluctuations observed during this uplift episode. We show that the conduit-dominated magma transfer between both reservoirs controls the dynamics of the system. A deep mafic reservoir recharging an evolved shallow reservoir would explain the mafic enclaves found in the dacites in the latter eruption and thus offers a physical model to jointly explain the observations we get from petrology, geochemistry and geophysics, bridging the disciplines. The model presented here can be used to study potential recharge mechanism occurring in the different stages of an eruptive cycle including pre-eruptive, co-eruptive and post-eruptive stages.

## 1. Introduction

Geodetic modeling inferred from InSAR and GNSS data has become the standard approach to interpret the mechanical processes governing

volcanic plumbing systems at depth. According to the classical volcanic deformation cycle, uplift episodes prior to eruptions are interpreted as magma migration to the shallow crust. In contrast, subsidence episodes are associated with the emptying of the reservoir during eruptions (Biggs

<footnote><sup>a</sup> COMET, School of Earth and Environment, University of Leeds, Leeds, UK</footnote>

<footnote><sup>b</sup> GET/UMR5563 (UPS, CNRS, IRD, CNES), Obs. Midi-Pyrénées, Université P. Sabatier, Toulouse, France</footnote>

<footnote><sup>c</sup> Centro Sismológico Nacional, Universidad de Chile, Santiago, Chile</footnote>

<footnote><sup>d</sup> Dirección de Desarrollo e Innovación, Universidad de Concepción, Victor Lamas 1290, Concepción, Chile</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: C.P.NovoaLizama@leeds.ac.uk (C. Novoa Lizama).</footnote>
and Pritchard, 2017). However, this interpretation does not explain the large and long-term uplift episodes observed at calderas that do not end in eruptions or volcanoes that erupt without prior uplift (Pritchard et al., 2019; Biggs et al., 2014; Wicks et al., 2011; Novoa et al., 2022). Additionally, there are volcanoes that exhibit no deformation during eruption and some of them even undergo episodes of uplift, which contradicts the classical model (Yip et al., 2022; Champenois et al., 2014). Here, we analyze the long-lasting uplift episode observed during the last eruption of the Nevados de Chillan Volcanic Complex (NdCVC) using long-term GNSS and InSAR data that span the eruptive period from 2015 to early 2022. Based on seismological and petrological studies at this volcano that evidence two magmatic chambers interacting in previous eruptions, in this study we explore the potential mechanism of recharge that would explain the uplift episode observed. We show how the interaction of these two chambers conditions the temporal evolution of these displacements and would explain their differences in the near and far field of the volcano. Such differences have been observed in some other volcanoes at different stages of their eruptive cycle, so the recharge mechanism proposed here is not limited to the study of observations only during eruptions, but during pre-eruptive or post-eruptive periods can also be applied.

The Nevados de Chillan Volcanic Complex (NdVCV) is a large composite stratovolcanic complex located in the Southern Andean Volcanic Zone, which extends roughly from 33° to 46° S and results from the subduction of the Nazca oceanic plate under the continental South America plate. This volcanic complex is an extensive 16-km long chain

of stratovolcanoes and several domes that are aligned along a NNW-trending axis (Fig. 1). This alignment, which is not expected from the current strain partitioning in the Southern Andean Volcanic Zone, has been proposed to be controlled by NW-SE trending left-lateral strike-slip fault inherited from pre-Andean geological processes (López-Escobar et al., 1995; Cembrano and Lara, 2009). The two main eruptive sub-complexes Cerro Blanco and las Termas are located to the north-west and the south-east ends of the NNW-SSE alignment and are distinguished by the different composition of their predominantly andesitic and dacitic eruptive products, respectively (Mee et al., 2009; Lupi et al., 2020; Oyarzun et al., 2022; Dunkley and Young, 2000; Dixon et al., 1999). NdCVC is one of most active volcanoes in Chile, ranked fourth by the Chilean Geological Survey in term of hazards (www.sernageomin.cl). Since the 17th century, volcanic activity has been mainly concentrated in the Las Termas subcomplex comprising four cones: Nuevo, Arrau, Viejo and Chillan, Fig. 1, and have evolved to a relatively uniform dacitic composition over the last century, associated to the relatively small eruptions occurred between 1906 and 2008 (Smithsonian Institution Global volcanism report available at http://www.volcano.si.edu; (Naranjo et al., 2004; Oyarzun et al., 2022; Coppola et al., 2016; Orozco et al., 2016).

The latest eruptive event at the NdCVC has been characterized by different volcanic activity linked to the hydrothermal system and the magmatic source. From 2012 to 2015, the activity was characterized by hydrothermal system predominance and low number of seismic events, without superficial activity until December 2015. On January 8, 2016 an

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Reference map of the study area in southern Chile correspond to the white box in the top right-hand inset map. The green circles and their associated name show the location of the five GNSS stations. The red polygons show the distribution of eruptive products during the Last Glacial Maximum period comprising activity in the periods 1906–1946, 1973–1986, 2008 and 2016-present. The main units are described as VN, Volcán Nuevo; VA, Volcán Arrau; VV, Volcán Viejo; VC, Volcán Chillan. Coordinates are expressed in UTM-WGS84 (19 zone South). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
eruption began with a strong explosion observed at the surface, followed by several explosive episodes occurring between January and August 2016. Predominant seismicity and the clear magmatic composition of the gas at the onset on the eruptive activity showed that the eruption was driven by pressurized magmatic gases, followed by eight months of volcanic activity linked to hydrothermal processes (Cardona et al., 2021; Moussallam et al., 2018). Then, a transitional phase with a domination of phreatomagmatic activity was identified until November 2017, which has evolved to a magmatic phase characterized by an increase of seismicity, including the presence of shallow LP events and the emission of four lava flows, which extruded a total eruptive volume of 300,000 m³ between August 2019 and December 2019. The first lava flow occurring in August 2019 was correlated with the onset of an uplift event starting at the end of July 2019 and which continued until the end of 2020 (Cardona et al., 2021; Astort et al., 2022). Since then, volcanic unrest has been persistent manifesting in intermittent effusive-explosive eruptive style, characterized by dome growth and destruction, small ash explosions, gas release and emission of pyroclastic and lava flows. After six years and 10 months of eruptive activity, on January 3, 2023, the Chilean volcano Observatory (OVDAS) has changed the alert to green, indicating that this volcano has returned to background activity and the eruption has ended (https://rnvv.sernageomin.cl/complejo volcanico-nevados-de-chillan/).

The presence of mafic enclaves in some dacites observed in the eruptive products of the 2016 eruption suggests the ascent of mafic magma from depth and mafic-felsic interactions in the shallow part of the plumbing system beneath NdCVC (Oyarzun et al., 2022; Moussallam et al., 2018). This interpretation is in agreement with the presence of a confined region of partial melt at depth beneath NdCVC revealed from ambient seismic noise Rayleigh-wave tomography (González-Vidal et al., 2018). The model derived from the tomographic inversion shows that the magmatic plumbing systems has a vertical extension of more than 15 km (González-Vidal et al., 2018). A recent petrological study by Oyarzun et al. (2022) using thermobarometry not only supports the results of the previous study but also advances our understanding of the structure of this underlying magma plumbing system. This study suggests a vertical compositional zoning in the feeder system spanning from 2 to 17 km of depth, which would explain the different compositions in the Cerro Blanco and Las Termas subcomplexes (Oyarzun et al., 2022). Therefore, denser magmas such as basaltic andesites and andesites originate from a deeper part of the magmatic system, at 7 to 17 km of depth. Conversely, evolved and less dense magma such as dacites and rhyodacites originate from shallower depths between 2 and 9 km. This zonation would be the result of magma upwelling from a MASH zone (melting, assimilation, storage and homogenization) located at the crust-mantle boundary and encountering various crustal discontinuities within the crust. Magma migration from the MASH towards the surface would explain the eruptive products from the Cerro Blanco subcomplex, as well as those in the Las Trancas subcomplex. This migration is also supported by the presence of basaltic andesitic enclaves in some dacites and rhyodacites from the last eruptions. It is from these observations, in conjunction with mineral chemistry and geochemical modeling, that it has been proposed that the recent eruptions in Las Trancas subcomplex were caused by mafic magmas coming from a deeper reservoir that re-activated a shallower evolved reservoir from which the dacitic material reached the surface (Oyarzun et al., 2022).

This study investigates the potential recharge mechanism associated to the crustal architecture of the NdCVC plumbing system evidenced by recent seismological and petrological studies, to understand the observed ground surface displacements captured by InSAR and GNSS from 2015 to early 2022 (Oyarzun et al., 2022; Moussallam et al., 2018, 2021). The feasibility of the two-chamber recharging process proposed for past eruptions during the last eruption of the NdCVC is assessed by analyzing the spatio-temporal behavior of this long-term displacement field. The locations and geometries of the two reservoirs are investigated using a combination of analytical and boundary element modeling,

whilst the associated recharge dynamics is studied using Bayesian inversion and the Ensemble Kalman Filter (EnKF) technique.

## 2. GNSS and InSAR time series

Ground deformation at NdCVC is characterized from daily solutions at five continuous GNSS stations installed by the Observatorio Volcanológico de los Andes del Sur (OVDAS) covering a period of six years from January 2016 to early 2022 (Fig. 2). The GNSS time series were decomposed into a linear combination of principal components by using the Principal Component Analysis-based Inversion Method (PCAIM) (Kositsky and Avouac, 2010) (see details in Appendix A.1.1). This analysis suggests that the first component of the decomposition explains 95 % of the total variance in the GNSS data. While this indicates that the displacements are mostly explained by the first component, there are systematic residuals that are significant. The first component of the decomposition fails to capture subtle variations in the horizontal component of ground surface displacements observed at various stations during the first six months of the acquisition period. If only the GNSS recordings acquired from the start of the uplift episode (July 2019) are considered, the first component of the decomposition well reflects the temporal evolution of the displacement recorded at the four stations located on the flank of the volcano, with a reduced chi squared of the residual, $\chi_r^2$, ranging from 0.22 to 0.29 and from 0.18 to 0.28 for the horizontal and the vertical components of the displacements, respectively. This result highlights the strong spatio-temporal correlation observed between the signals recorded by these stations. During the same uplift period, the first component does not fully capture the temporal pattern of the displacement observed at station TRAN (Fig. 2) as evidenced by $\chi_r^2$ values of 0.8 and 0.85 for the horizontal and vertical components of the displacement, respectively. The main conclusion of this analysis is that even it is theoretically possible to mimic the observed displacements within the uncertainty using a single-source model, it will be impossible to fit the subtle variation of the displacement observed at TRAN station.

SAR imagery at NdCVC was acquired by two satellite missions: C-band Interferometric Wide Swath mode images from the European Space Agency (ESA) Sentinel-1 A/B and L-band Stripmap mode images from ALOS-2 operated by the Japan Aerospace Exploration Agency (JAXA). All treatment details of the InSAR images are detailed in Appendix A.1.2. Fig. 3 shows some wrapped InSAR interferograms spanning different time intervals during the 2015–2022 period.

We assessed the consistency between GNSS and InSAR at geodetic monitoring sites in the study area. Both times series are in good agreement with the LOS-projected GNSS measurements with a rms of 0.4 and 0.7 cm for ALOS-2 and Sentinel-1 data, respectively (Fig. 3 and Appendix A.1.3).

## 3. GNSS and InSAR time series analysis during the eruption

Time series of GNSS and InSAR (Fig. 2–4) highlight the wide range of temporal behavior of ground surface displacements observed from 2016 to 2022 at NdCVC during the eruptive period. In this section we compare the temporal evolution of these ground surface displacements at NdCVC with the volcanic activity information collected from different studies and from the public reports generated by the Volcanological observatory OVDAS (https://rnvv.sernageomin.cl/complejo-volcanico-nevados-de-chillan/). A schematic diagram of the temporal relation between volcanic activity and the vertical displacement observed at TRAN and FRSC is shown in Fig. 4.

### 3.1. From no deformation to slight subsidence between January 2016 to June 2019

Until August 2016, volcanic activity was mainly characterized by phreatic eruptions driven by magmatic gases with several explosive
<figure><img src="image_6.png" /><figcaption>Fig. 2. GNSS time series from January 2016 to December 2021. The blue points are the GNSS measurements after reducing the post-seismic and seasonal displacement using the function fitting models derived from GNSS measurements recorded at ANTC. Red lines show the reconstructed GNSS displacements using the first principal component of the PCA decomposition. The data variance explained by the first component is 95 % for each one of the components of the GNSS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

episodes occurred between January to August 2016 (Cardona et al., 2021; Global Volcanism Program, 2023). During this period no vertical displacements are observed (P1 in Fig. 4 and Fig. 2) but, as mentioned above, subtle variations in the horizontal component of ground surface displacement are visible. As these variations are not associated with vertical displacement, the decomposition using the first principal component fails to capture them (Fig. 2). In the Sentinel-1 interferogram spanning the period March 2015 to March 2016, a subtle interferometric signal is also observed but we cannot rule out this being related to tropospheric effects rather than ground surface displacements (Fig. A.11).

The phreatic period was followed by a series of phreatomagmatic eruptions, accompanied by high energy explosive eruptions with a moderate presence of juvenile components and the emergence of a new crater (named Nicanor) to the east of the volcano Nuevo in October 2017 (Lupi et al., 2020; Cardona et al., 2021), P2 in Fig. 4. From December 2017 to November 2018, a lava dome grew by the extrusion of viscous magma from a volcanic conduit. This lava dome grew gradually, aggregating a maximum volume of about 370,000 m³, which was partially destroyed in August 2019 (Cardona et al., 2021). GNSS time series spanning the phreatic period, from October 2017 to August 2019, show a slight subsidence recorded by all GNSS stations with the single exception of TRAN (P3 in Fig. 4). This slight subsidence prior to the uplift period also seems to be observed in the InSAR Sentinel-1 time series, although the precision of the InSAR measurements makes it impossible to state it with certainty (Fig. 3, Fig. A.12). A maximum subsidence of 0.8 cm y⁻¹ is observed at the FRSC station between the start of the vertical displacement time series in January 2016 and before the onset of the uplift period in July 2019. One possible explanation for this subsidence could be that our correction of the post-seismic

contribution may be slightly excessive. However, if this were the case, this subsidence would be observed in all the stations to which we applied the correction, but it is only observed at stations located on the volcanic edifice and not at TRAN station, located to the southwest of the volcanic complex (Fig. 1). This leads us to conclude that it is more likely to be caused by a process related to volcanic activity, such as a decrease in pressure in a superficial hydrothermal system or in a shallow magma reservoir located underneath the volcanic edifice.

## 3.2. Uplift episode observed between July 2019 to January 2022

Following the partial destruction of the lava dome, a strong effusive activity period of around four months began at the end of August 2019 with four pulse of lava flows descending down the northeast flank accumulating a total of 300,000 m³ of extruded material (Cardona et al., 2021). The first lava flow occurring in August 2019 was correlated with a significant increase in seismic activity (most of them LP earthquakes) and explosions per day, as well as with the onset of an inflating event starting at the end of July and observed over the 5 GNSS stations installed at NdCVC, as well as in InSAR observations (P4 in Fig. 4, Fig. 2–3). One of the most striking aspects of this uplift event is the fast increase of ground displacements at the beginning, followed by a slow decay over the next years (Fig. 2–4). This temporal evolution of the displacements has been observed at different volcanoes and has been characterized by a first order exponential decay (Rodríguez-Molina et al., 2021). Similarly, the first component of the decomposition of GNSS displacements observed during this uplift event at NdCVC can be approximated using a decreasing exponential function as follows:

$$a(1 - e^{-t/\tau}) + b, \qquad (1)$$
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Mosaic of a subset of ALOS-2 and Sentinel-1 interferograms and time series of GNSS and InSAR data at FRSC. Left) The satellite to ground radar line of sight (LOS) is shown with black arrows. A complete cycle of phase (blue, yellow, red) represents a decrease in range of 2.8 and 11 cm between the ground surface and the satellite for Sentinel-1 and ALOS, respectively. Areas without coherence are shown in gray. Right) Red points are the LOS-projected GNSS data. Blue points are the InSAR data with their standard errors. Due to the lack of coherence no Sentinel-1 observations are available at the location of the GNSS station TRAN. Rms between the LOS-projected GNSS measurements and ALOS-2 and Sentinel-1 data are 0.4 and 0.7 cm, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Schematic diagram of the temporal relation between volcanic activity and the vertical displacement. Red and blue dots show the vertical displacements observed at TRAN and FRSC, respectively. Black line represents the number of LP per day extracted from Cardona et al. (2021). P1 Predominance of Phreatic activity driven by magmatic gases; P2 Transitional stage ash emission occurrence of explosion; P3 Lava Dome building and the partial destruction of the dome and its transitory constructions; P4 New magma input and lava flows; P5 Explosions and pyroclastic flows continue new dome emerges from Nicanor crater Dome construction (15/11/2020); P6 Pyroclastic flows, gas-and-ash plumes, and a new lava dome during December 2021-May 2022; P7 about 1.5 million cubic meters of effusive material was emitted during 17 September-2 December 2021. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
where the coefficient ‘a’ characterizes the amplitude, ‘b’ is an offset shifting the curve vertically and ‘r’ is a time constant. In the case where all the displacements recorded by the GNSS stations can be explained by a single component, this indicates that only the coefficients a and b vary depending of the location of the GNSS stations but tau is constant whatever the GNSS station considered. This is because the ground displacement field is spatially stationary and that only its amplitude varies over time. This is indeed what we observe for the four stations located on the flank of the volcano characterized by the same exponential decay constant of 0.76 yr with a 95 % confidence interval of 0.74–0.78 yr (Fig. 5.B). The only station that cannot be explained by this decay constant is the TRAN station, which is located southwest of the volcanic complex (Fig. 5.B). Indeed, at TRAN station the best-fit value for the exponential deformation trend is about twice that observed on the other station with a value for the exponential decay constant of 1.51 yr (1.41–1.62 yr at 95 % confidence).

The main result to be drawn from the analysis of the uplift episode observed during the magmatic phase of the eruption is that the four stations located on the volcanic edifice are experiencing the same pattern of temporal evolution, but with different amplitude depending on the location of the stations. This leads us to consider that ground displacement observed at these stations is related to the same source located beneath the volcanic edifice, which would be pressurising during this period. On the other hand, the different temporal behavior of displacements recorded by the TRAN station suggests that another source is contributing to its observed ground displacement. From these observations and in line with petrological and geochemical studies, we hypothesize that the surface displacements observed over the volcanic complex are related to the dynamics of two sources. In such a context,

changes in pressure in a shallow reservoir would be the main contributor to the displacements observed in the study area. However, TRAN would have also recorded displacements triggered by pressure changes in a deeper reservoir that is recharging the shallow one.

## 4. The dynamics on a double-reservoir plumbing system: modeling approach

To constrain the geometry, location and potential interaction between the hypothesized two reservoirs in this volcanic complex, this section describes the double-reservoir model used in this study. The aim is to replicate the temporal behavior of surface displacements observed at NdCVC during the observed uplift episode of the magmatic phase, spanning from June 2019 to January 2022.

### 4.1. Analytical solution for a double-reservoir model accounting for any reservoir geometry and topography

Reverso et al. (2014) proposed an analytical model describing the dynamics between two magma reservoirs at different depths connected by an incompressible-magma-filled hydraulic pipe. This model considers two possible geometries for the reservoir: a spherical or a sill-like source, which are embedded in an homogeneous elastic half-space, with a flat free surface. In this study, we derive an analytical solution based on a similar mass balance of the system, but our model allows for flexibility in the geometry of both reservoirs and accounts for the topography in the zone. We develop this model considering that our study area exhibits a significant relief (altitude ranging from 800 to 3200 m.a.s.l), preventing the use of classical half-space models (Cayol and Cornet, 1998). A conceptual diagram describing the double reservoir model adopted here

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Time function of the first component and vertical displacement time series spanning from July 2019 to January 2022. A) Time function of the first component from the decomposition. Red line shows the resulting best fitting exponential function with an exponential decay constant of 0.76 yr (0.74–0.78 yr). B) The black points represent the vertical displacements observed at the GNSS stations. The exponential decay constant is too high for a first-order fit to the displacements observed at TRAN, unlike at other stations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
is shown in Fig. 6.

In this model the evolution of pressure changes in the shallow ($\Delta P_s$) and the deeper reservoir ($\Delta P_d$) are controlled by the geometry of both reservoirs, plus six parameters: magma viscosity ($\eta$), magma density ($\rho_m$), ($\rho_c$) is the difference of density between the magma and the bedrock i.e. ($\rho_r-\rho_m$), radius of the cylindrical conduit ($a_c$), $H_c$ is the difference of depth between the deep and shallow reservoirs, i.e. $H_c = H_d - H_s$, and a constant magmatic flux ($Q_{in}$). Note that in this model we assume that the uplift observed during this period is due to material arriving from a deep chamber recharging a shallow chamber. Given the absence of subsidence and the lack of an estimate of the eruptive volume available during the period June 2019 to January 2022, we do not consider a magmatic flow to the surface in our model. However, we will discuss how this assumption might affect our calculations in section 5.2 below.

From the mass balance detailed in Appendix A.2, we obtain the following first-order linear equation that describes the dynamic of the system:

$$ \frac{d\Delta P_s(t)}{dt} + \xi \Delta P_s(t) = \frac{\pi a_c^4}{C_s 8 \eta H_c} \left[ gH_c \rho_c + \frac{Q_{in} t}{C_d} + \frac{C_s}{C_d} \Delta P_s(0) + \Delta P_d(0) \right] \quad (2) $$

with

$$ \xi = \frac{\pi a_c^4}{C_s 8 \eta H_c} \left( 1 + \frac{C_s}{C_d} \right) $$

where $\Delta P_s(0)$ and $\Delta P_d(0)$ corresponds to the initial overpressure at the shallow and deeper reservoir respectively and $C_s$ and $C_d$ are constants that describe the relationship between pressure change and volume

change in both reservoirs as following:

$$ \begin{aligned} \Delta V_s(t) &= C_s \Delta P_s(t) \\ \Delta V_d(t) &= C_d \Delta P_d(t) \end{aligned} \quad (3) $$

The solution of this first-order linear differential equation for the pressure evolution in the shallow and deeper reservoirs can be then expressed as:

$$ \begin{aligned} \Delta P_s(t) &= \frac{Q_{in}}{C_d + C_s} t + \Delta P_s(0) + A(1 - e^{-\xi t}) \\ \Delta P_d(t) &= \frac{Q_{in}}{C_d + C_s} t - \frac{A}{C_d}(1 - e^{-\xi t}) + \Delta P_d(0) \end{aligned} \quad (4) $$

where,

$$ A = \frac{C_d}{C_d + C_s} \left( \rho_c g H_c - \Delta P_s(0) + \Delta P_d(0) - \frac{8 \eta H_c C_s Q_{in}}{\pi a_c^4 (C_d + C_s)} \right) $$

The constants $C_s$ and $C_d$ depend on the shear modulus and the geometry of each reservoir. We express the linear relationship between surface displacement and the pressure variations in the chambers of the dynamical model using a matrix **H**.

## 4.2. Inversion strategy

To perform the inverse modeling of the volcano deformation, we use Defvolc, a 3D Boundary Element Method (BEM) for elastic media combined with a non-linear inversion method. Defvolc is a freely available software accessible online (http://www.opgc.fr/defvolc/Vue/MainPag e.php). The numerical modeling method can incorporate realistic

<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Left: Representative model of the architecture deducted below NdCVC from the most recent petrological study in the area (Oyarzun et al., 2022). The vertical compositional zoning is shown by the yellow area representing dacites/rhyodacites and the blue area indicating the depth of basaltic andesites-andesites. Center: Conceptual diagram of the model described in section 4. Right: Schematic steps indicating the modeling strategy used in this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
topography and any number of pressure sources with different geome-
tries. Throughout all the stages of our modeling, we assume that the
reservoirs are embedded in an isotropic homogeneous elastic half space
limited by the surface topography and that a uniform pressure is exerted
over the whole boundary of the sources. For all the models, we assume a
Poisson coefficient (μ) of 0.25 and a shear modulus (G) of 4 GPa (Heap
et al., 2020).

The inversion strategy we adopted to search for the geometry of the reservoirs is as follows: using cumulative displacements of GNSS and InSAR data. In Defvolc, we first estimate the parameters of the shallow source, whose contribution to ground displacement is dominant. We then subtract the contribution from the shallow source to find the deep source that best fits the remaining observed displacements. This approach is based on the result of the principal component analysis, which shows that the ground displacements can be replicated within uncertainties using a single-source model. The identification of the optimal geometry for the shallow source allows to find the value of C<sub>s</sub>, from the values of ΔP<sub>s</sub> and ΔV<sub>s</sub> given by Defvolc. Then, with the geometry found, the surface displacements caused by an over-pressurisation of 1 MPa are simulated over the source, which allows us to calculate the value of **H** for the shallow source. Next, for the deep source we proceed similarly, calculating the C<sub>d</sub> constant and the value of **H** for the deep source, both of which are derived from the Defvolc inversion.

4.3. Shallow magma chamber

Multiple lines of evidence suggest that the shallow source beneath the volcanic complex is located at depth between 3.5 and 7 km (Cardona et al., 2021; Oyarzun et al., 2022; Astort et al., 2022). In agreement with that, we therefore consider a variety of classical source geometries (sphere, spheroid, horizontal and circular crack) to model the surface displacement observed by GNSS and InSAR constraining the inversion to search for sources in such a depth range. We performed a joint inversion of the cumulative displacements from January 2019 to May 2020 using the ascending ALOS-2, descending Sentinel-1 and GNSS data. Inversions favor a prolate spheroid located at 5.8 km depth below the Earth surface elongated in the direction N126°E with a major axis of about 3.3 km, slightly dipping to the East by 12°. The volume of the source is about 2.84 km³ and the estimated geodetic volume change of the source is about 27.5 ±6 10⁶m³. This best model yields χ²_r of 9.2 and 3.2 for the horizontal and vertical components of GNSS observations, respectively. For InSAR data χ²_r are 1.20 and 0.59 for descending Sentinel-1 and ascending ALOS-2, respectively. The model reproduces correctly the global displacement pattern but falls short in fully replicating the displacements observed by both InSAR and GNSS measurements in the summit region of the volcanic complex (Fig. 7). This explains in particular the χ²_r of 9.2 obtained for the horizontal components of the GNSS. Despite this localized high residual, the prolate spheroid source embedded in an elastic medium seems a reasonable approximation for modeling the surface displacements, although the shape of the source and the rheology of the surrounding crust are certainly more complex than those assumed in our model. In the next section we explore whether the inclusion of a deep magmatic chamber connected to this shallow chamber can improve the fit to the data either spatially or temporally.

4.4. Deep magma chamber

As mentioned in section 3.2, we suspect that the peculiar behavior of
the displacements observed at the TRAN station could be due to varia-
tions in the internal pressure of a deep magma reservoir that is
recharging the shallow one. Therefore, in agreement with petrological
and geochemical observations (Oyarzun et al., 2022), here we search for
the overpressurization on a deep source beneath the southwest part of
the volcanic complex able to explain these observations. The significant

residual between the observed and modeled displacements using the
prolate spheroid, particularly for the horizontal component of the
displacement, and the lack of surface displacements observed with
InSAR in the far-field, prevents us from being able to estimate the deep
source geometry from the signal that remains after subtracting the
contribution of the shallow source. Therefore, we assume the deep
reservoir to be a planar circular source, similar to a sill-like source. Then
the double-reservoir model described in section 4.1 is used to assess the
temporal evolution of its pressure to explain the temporal evolution of
the displacements registered in the GNSS stations.

4.5. Inversion to constrain the double-reservoir dynamics

To reproduce the temporal evolution of displacements at NdCVC, we
only consider the vertical displacements acquired by GNSS at TRAN and
FRSC. This choice is based on the following points: First, the TRAN
station potentially reflects the highest contribution of the deep source at
the surface. Secondly, the discrepancy (χ²_r of 9.2) between the observed
horizontal displacement components and those predicted by the shallow
source model prevents their use. Additionally, due to the strong spatio-
temporal correlation of the vertical displacements recorded by the sta-
tions located near the top of the volcano, we chose FRSC as it has the
strongest observation potentially related to the shallow source.

The vertical components of GNSS are inverted using a Bayesian-
based inversion through the Markov Chain Monte Carlo (MCMC) algo-
rithm using the double-reservoir model described in section 4.1. The
geometry and the depth of the shallow source, as well as the values of C_s
and the vertical displacements used in the inversion are computed from
the 3-D mixed BE method DefVolc. First, we considered 6 model pa-
rameters to be determined by GNSS data inversion (Table A.2). The
minimum and maximum bounds the magma viscosity and the difference
of density between the magma and the bedrock (ρ_c) are chosen by
considering several studies (Takeuchi, 2011). The sensitivity analysis
presented in Appendix A.3 indicates that there are several strong trade-
offs between different parameters that do not enable the inversion
process to determine simultaneously these parameters. To eliminate the
trade-off between magma viscosity (η) and the radius of the conduit (a_c)
we fixed η to a mean value of 10⁴ Pa·s. In the same way to eliminate the
trade-off between the depth of the deep source and the difference of
density between the magma and the bedrock we fixed the value of ρ_c to
2500 kg·m⁻³. The resulting 1-D and 2-D marginal posterior probability
density functions deduced by the GNSS data inversion are plotted in
Fig. A.14.

4.5.1. Constant basal flow (Qin)

First, we consider all the model parameters constant over time.
Table 1 summarize the lower and upper bound search value for the
uncertain model parameters and reports the best model parameters with
uncertainty level of 95 % deduced by the inversion. The resulting two-
magma reservoir model inferred from GNSS data matches well the
observed 900 days-long time series of ground displacements (Fig. 8A).
The best value of the radius conduit a_c is 1.82 m (CI at 95 % of 1.79–1.85
m), and the one of the constant basal magma flow (Q_in) is 0.0156 km³y⁻¹
(CI at 95 % of 0.0140–0.0171 km³.y⁻¹). During this period the shallow
reservoir experiences a volume increase of 26x10⁶m³ with a pressure
increase of 38 MPa (Fig. 4A). Such a volume change accounts for about
97 % of the uplift episode observed during the 900-days period. On the
other hand, the volume change in the deep reservoir decreases to
–5.8x10⁶m³ the first 230 days, and then increases until 12x10⁶m³
(Fig. 8A). The decrease in volume change in the deeper reservoir during
the first 230 days then indicates that the material flowing through the
conduit into the shallow chamber during this period exceeded the
recharge associated with the basal magma flux (Q_in) of 0.0156 km³y⁻¹.
In contrast, after these 230 days, the volume increase in the deep
reservoir indicates that material flowing into the shallow reservoir from
<figure><img src="image_12.png" /><figcaption>Fig. 7. Modeling result for the prolate spheroid source. Top row) ascending ALOS orbit; LOS displacements towards the satellite are negative. Middle row) Descending Sentinel-1 orbit; Bottom row) GNSS data.</figcaption></figure>
<table><caption>**Table 1**
Parameters of the two-magma model inferred from the inversion of GNSS data.</caption><thead><tr><th>Parameters</th><th>Min. Values</th><th>Max. Values</th><th>Lower (95 %)</th><th>Upper (95 %)</th><th>Best Values</th></tr></thead><tbody><tr><td>a<sub>d</sub> (km)</td><td>0.1</td><td>5</td><td>1.50</td><td>2.49</td><td>1.84</td></tr><tr><td>a<sub>c</sub> (m)</td><td>0.5</td><td>10</td><td>1.79</td><td>1.85</td><td>1.82</td></tr><tr><td>H<sub>d</sub> (km)</td><td>8</td><td>30</td><td>14.50</td><td>15.43</td><td>15.94</td></tr><tr><td>Q<sub>in</sub>(km3.y<sup>-1</sup>)</td><td>10–3</td><td>10–1</td><td>0.0140</td><td>0.0171</td><td>0.0156</td></tr><tr><td>ρ<sub>c</sub> (kg.m<sup>-3</sup>)<sup>*</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>250</td></tr><tr><td>η (pa.s)<sup>*</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>10<sup>4</sup></td></tr></tbody></table>

* Indicates the value of the parameter is fixed in agreement with the conclusion of the sensibility analysis presented in Appendix A.3.

the deep reservoir decreased, rising at a lower rate than the basal magma flux ($Q_{in}$) of 0.0156 km³y⁻¹ that was recharging this deep chamber. This slower rate of material transport would then explain the slower rate of volume increase in the shallower chamber after the initial 230 days and therefore the slower rate observed in the displacements at FRSC (Fig. 8A).

Based on these results, our model indicates that the observed uplift at NdCVC can be explained by constant basal magma flux ($Q_{in}$) of 0.0156 km³y⁻¹. However, it is the variation of magma transport in the conduit connecting the two chambers that conditions the dynamics in the system. This variation of transport rate from faster to slower is consistent with a first exsolution of gases opening the conduit driving magma towards the upper chamber at a high pressure and then giving way to a degassed magma rising more slowly.

To assess whether this deep reservoir improves the spatial fit of the data, we include this double reservoir model to fit the Sentinel-1 and ALOS-2 data in Fig. A.13. We can see that by including this deep reservoir, we have a larger deformation area, affecting the TRAN station, however, the magnitude of the displacements produced by the deep reservoir is one order lower (~1 cm) than that produced by the shallow chamber (~10 cm). Furthermore, we would like to note that with the double chamber model used here, we still have a strong residual in the summit region of the volcanic complex (Fig.A.13).

### 4.5.2. Variable basal flow ($Q_{in}(t)$)

Although our model captures well the long time behavior associated with an exponential decay over the three years of the magmatic phase of the eruption, it cannot reproduce the small variations observed in FSRC, Fig. 8A. Here we test whether these temporal variations are related to change of the magma supply rate at depth and therefore the assumption of constant basal flow rate ($Q_{in}$) in our model in section 4.5.1 would be a strong limitation to reproduce these variations. To test this hypothesis, we perform an Ensemble Kalman Filter (EnKF) approach similar to the one proposed by Bato et al. (2018), which description in details is given in Appendix A.4. In this approach, the governing eq. A.24 is used to find the best variation of $Q_{in}$ that explain the variations in the vertical displacements observed between the 2019 and 2022.

The EnKF results show that variations of magma supply rate around its mean value of 0.0156 km³y⁻¹ well explains the subtle variation of vertical displacements observed at the GNSS stations located on the summit of the volcano and at TRAN station (Fig. 8B). Similar to the model in section 4.5.1, here we show a decrease in volume during the first 230 days in the deeper reservoir, indicating a transport of magma from the deep to the shallow reservoir higher than the basal magma flow into the deeper reservoir from the MASH ($Q_{in}$). Nevertheless, a difference in behavior is observed after these 230 days, particularly in the period covering the 550 days to the 670 days. We show that the deeper reservoir pressure is not increasing for all the remaining period as in the previous model (Fig. 8A), instead it decreases at 550 days and then increase again at 670 days (Fig. 8B). This slight decrease in pressure would be associated to a lower magma basal flow into the deeper reservoir ($Q_{in}$), therefore, the volume of material transported into the shallow

chamber would be greater than the material entering the deep reservoir ($Q_{in}$) in this period. These slight differences in the magma inflow ($Q_{in}$) would then explain the small uplift accelerations observed in the different GNSS stations at the summit of the volcano, Fig. 8B.

## 5. Discussion

### 5.1. Phreatic and phreatomagmatic phases of the eruption explained by hydrothermal activity and degasification of the shallow reservoir

The onset of the 2016 eruption was evidenced by the first explosion observed at the surface of the NdCVC in January 2016. The characteristics and evolution of the volcanic activity during the first eight months suggest that the volcanic emission was mainly linked to hydrothermal processes. However, several lines of evidence suggested that this phreatic activity was mainly driven by a magma intrusion that lasted only a month, as a high percentage of juveniles components was evident only until February 2016, and the long period events were located at less than 5 km below NdCVC that first month (Cardona et al., 2021; Moussallam et al., 2018; Benet et al., 2021). According to these observations, Moussallam et al. (2018) proposed that the eruption was caused by the recharge of a shallow reservoir fed by a deep source, suggesting that a similar mechanism could have been occurring during the past eruptions at NdCVC. Although all these observations point to a short-period magmatic intrusion, no episode of uplift was detected either by GNSS or INSAR during this phreatic phase. Continuous GNSS observations recorded at the stations located on the volcanic edifice only reveal subtle horizontal displacements (Fig. 2), whose origin could be associated to fracture opening. This lack of an uplift episode suggests that the volume change induced by the magma ascent during this magma intrusion was too small to engender detectable ground displacement, but enough to reheat the hydrothermal system. On the other hand, the addition of gases, heat and/or mass provided by this intrusion's mafic magma would promote the remobilisation of the viscous magma stored in the shallow reservoir upon reaching it. This remobilisation inside the shallow chamber would then explain the slight subsidence highlighted by GNSS and INSAR observations from September 2016 to August 2019, which roughly correlates with the extrusion of a lava dome, with a total volume of about 370,000 m³, and its destruction. The extremely slow growth of the dome observed in the NdCVC (Moussallam et al., 2021) would then not be associated with a magma recharge in this period, but could rather be associated with the slow transport of this small amount of magma from the intrusion to the surface, where the high viscosity of the magma in this shallow reservoir would prevents its fast ascent. Therefore, we suggest that the lack of deformation followed by the slight deflation observed previous to June 2019 would be related to the hydrothermal system, the degassing processes triggered by this short-period magmatic intrusion and its ascent towards the surface, Fig. 9. This relationship between the lack of deformation and degassing processes during eruptions has previously been observed at other volcanoes, where it has been explained by the high compressibility of a magma due to the large amount of volatiles it contains (Yip et al., 2022; Wasser et al., 2021). In terms of the seismic activity recorded over this period, we consequently suggest that this is mainly triggered by exsolved volatile flux rather than the magma flux into the shallow reservoir. Such a pattern of strong seismicity associated with no deformation or deflation during a dome extrusion event is not rare, as it has been observed on other volcanoes, in particular during the eruption of Soufrière Hills volcano at Montserrat (Mattioli et al., 1998; Christopher et al., 2015).

### 5.2. The dynamics of the double reservoir system during the magmatic phase of the eruption

The two-reservoir dynamical model we have developed for this study can well reproduce the spatial and temporal evolution of the uplift episode observed from July 2019 onwards. In line with previous studies
A) Constant basal flow

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

B) Variable basal flow

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

**Fig. 8.** Volume changes in the shallow and deep reservoirs are compared under constant and variable basal flow (Qin), together with the observed and modeled vertical displacements at the FRSC and TRAN stations. The red line represents the mean value and the gray fill is the spread (1σ). The red dashed lines indicate the contribution of the volume change of the deep source to the ground displacement captured by the GNNS stations. A) Constant basal flow, from the inversion of GNSS data using a MCMC algorithm. B) Variable basal flow, from the inversion of GNSS data using an Ensemble Kalman Filter approach. The yellow rectangle highlights the difference between both models: the deep reservoir pressure decreases after 550 days and then increases again at 670 days, which correlates with a lower magma basal flow into the deeper reservoir (Qin) during this period. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

that conclude there was a deeper reservoir feeding a shallow one during
the past eruptions, our model consists of an elongated spheroid reservoir
at 5.8 km depth below the surface connected to a deeper sill-like one
located at 15.4 km depth linked by a hydraulic pipe with a radius about
2 m. The largest contribution to ground surface displacement is

obviously related to pressure variation in the shallow reservoir, whereas
the deeper reservoir contributes to the surface displacements observed
far from the volcanic edifice where the shallow reservoir is located.
From the inversion of InSAR and GNSS data, this elongated spheroid
shallow source agrees with the geometry and location inferred from
<figure><img src="image_15.png" /><figcaption>Fig. 9. Conceptual scheme about the different stages observed during the eruptive cycle at NdCVC.</figcaption></figure>

surface displacements observed during the first year by Astort et al. (2022). The source is horizontal and its orientation is similar to the volcano chain with a strike angle of N126°E. This orientation suggests that its formation could be controlled by the same inherited NW-SE trending left-lateral strike-slip fault systems as the volcanic chain located above. The shallow depth of this shallow source makes it impossible to rule out that the origin of the uplift episode could be related to an increase of volume in the hydrothermal system beneath the volcanic complex. Nevertheless, considering that the uplift episodes were correlated with lava flows observed at the surface during the magmatic phase of the eruption and that the hydrothermal system in this zone is located further to the south east of the deformation source that we observe (Farias et al., 2014), our results support the hypothesis that a magmatic reservoir beneath NdCVC is the origin of this deformation.

We estimate the volume change of this shallow reservoir was 26 x10⁶m³, whereas this deep source ranges from -5.8 to 12x10⁶m³. These small changes in volume at depth explain why the surface displacements related to the deep source (less than 1 cm at TRAN) were not detectable in InSAR, see Fig. 8 and Fig. A.13. The decrease in volume in the deeper reservoir reaching 5.8x10⁶m³ during the first 230 days, would be explained by a magma flow from the deep to the shallow source that is larger than the volume entering the deep chamber from the MASH (Qin), Fig. 9. On the other hand, the increase in volume for the remainder of the time reaching 12x10⁶m³ indicates that during this second period, the volume of magma transported from the deeper reservoir to the shallow one is at a slower rate than the magma flow entering the deep chamber from the MASH (Qin) (Fig. 9). Astort et al. (2022) has reported a volume change of 44 +/- 14 x10⁶m³ between June 2019 and November 2020, associated to a shallow source that correlates with the shallow reservoir of our system. Our results suggest a lower volume change with a volume increase of 38 +/- 10x10⁶m³ between June 2019 and early 2022 associated with both sources.

Our estimate of volume change is limited by the model we used, which does not account for material flowing from the shallow reservoir to the surface. As indicated in Section 4.1, we did not include a conduit connecting the shallow reservoir to the surface in our model because there was no way to constrain it in the absence of observed subsidence during the period June 2019 to January 2022. Furthermore, the lack of a time series of volume fluxes defined in terms of dense rock equivalents spanning the eruptive period prevents the estimation of a plausible ratio of eruptive volume to geodetic volume change at NdCVC. The only

available estimate is the one reported by OVDAS of a volume of 12x10⁶m³ of erupted volcanic products until 2022 due to lava flows and domes, excluding the material erupted during the explosive activity (Pedreros-Delgado et al., 2024). Taking this value of 12x10⁶m³ into account, the estimated volume of magma flowing from the MASH of 38x10⁶m³ estimated in this study could be underestimated by about 30 %. Further studies estimating the total eruptive products, along with the percentage of gases and petrological data are needed to better constrain this volume change in the future.

In terms of the temporal evolution of ground surface displacements, our model assuming a constant basal flow (Qin) well reproduces the long temporal behavior associated to an exponential decay over the three years of the magmatic phase of the eruption. Nevertheless, it cannot reproduce the small variations observed at FSRC, Fig. 8A. As mentioned previously, the temporal evolution of the displacements observed by all the stations near the summit of the volcano are very similar. One possible explanation for these variations over the last three years of the eruption could be related to the rise of magma through an open conduit without producing large-scale surface displacement, and therefore explaining why these variations are not captured in the TRAN station. Nevertheless, magma ascent would have occurred without altering significantly the pressure in the deeper reservoir, since we did not observe these variations in the surface displacements recorded by TRAN. This last assumption seems inconsistent with our hypothesis of an interconnected magmatic system. Another possibility would be that these temporal variations are related to change of the magma supply rate at depth.

Accounting for an Ensemble Kalman Filter (EnKF) approach to assess variable basal flow (Qin(t)) we find that variations of the magma supply rate around its mean value of 0.0156 km³·y¹ can explain the subtle variation of vertical displacements observed at the GNSS stations located on the summit of the volcano and at TRAN station (Fig. 8B). In comparison with the constant basal flow model, an interesting difference of this model can be seen between the 550 and 670 days after the uplift started, where the decrease in magma inflow (Qin) would result in an increase in the transport of material to the shallow reservoir in this period explaining the small uplift accelerations observed in the different GNSS stations at the summit of the volcano, Fig. 8B.

Although our model well reproduces the temporal evolution of the displacements at NdCVC during the magmatic phase of the eruption, we cannot rule out the possibility that other mechanisms may also explain
the exponential decay of the displacement rate observed from 2019 to 2022. For example, some authors have considered the role of visco-elastic relaxation to explain the exponential decay observed in several volcanoes before, during or after eruptions (Dragoni and Magnanensi, 1989; Novoa et al., 2019; Del Negro et al., 2009; Newman et al., 2006)). On the other hand, Got et al. (2013) demonstrated that this temporal behavior of the displacement could also be explained by the effect of increasing damage around a magma chamber. Similarly, Segall (2016), has also shown that such a decay in an uplift event, could be the response of a sudden decompression in a magma of high compressibility contained in a reservoir of low compressibility. Although these different models could fit the temporal evolution observed at NdCVC equally well, the model proposed here best agrees with the geophysical and geochemical observations, which evidence a mechanism of recharge involving the contribution of two reservoirs to last eruptions at NdCVC. Furthermore, this model offers a simple way to characterise the potential recharge process acting in the volcanoes, allowing for the simulation of any type of spatial pattern observed in surface displacements, while also accounting for the topography in the models. Despite here we use this model to study the recharge mechanism occurring during the eruption, this model can be applied to study uplift episodes in other volcanoes that do not necessarily ends in eruptions. Finally, we think that this work highlight the importance of a detail analysis of the residuals between the models we make and the observations, since probably part of them could be associated with second or third sources dynamically acting during an eruptive event.

## 6. Conclusion

Based on petrology and seismology studies evidencing the existence of two interacting reservoirs beneath the NdCVC, in this study we elucidate the mechanism of recharge occurring during its last eruption from the analysis of a combined 6-year time series of InSAR and GNSS data. The model proposed couples an analytical approach with an open-source boundary element method to quantify the dynamics between these two reservoirs, taking into account the steep relief of the study area.

From the analysis of the surface displacements, our results show no deformation and a slight subsidence during the first three years of the eruption, which was followed by an uplift episode extending from June 2019 to January 2022. We propose that the first three years of no deformation and slight subsidence are associated with activity in the hydrothermal system and magma remobilisation in a shallow reservoir, caused by a short-term magma intrusion. We then propose that a new and larger magma intrusion occurs in June 2019, which continues for the next three years of the eruption decaying exponentially. Our model indicates that this intrusion was triggered by magma flowing from the MASH into the deep reservoir at a constant rate of 0.0156 km³/yr, with small changes in this constant rate explaining the small fluctuations observed during this episode of uplift. The uplift episode observed can be 97 % explained by the change in volume of a shallow reservoir at 5.8 km of depth, whereas the change in volume of a deeper chamber at 15 km depth accounts for the remaining 3 %. The shallow source dominates the temporal evolution of the four GNSS stations located in the volcanic edifice, whereas the deeper reservoir would explain the different temporal evolution of the surface displacements observed in the GNSS station located at 10 km of the volcanic summit. These displacements triggered by volume changes in the deeper reservoir were not detected by InSAR because they contribute less than 1 cm.y⁻¹ and were only detected by the continuous GNSS network at NdCVC. Therefore, this GNSS network was indispensable to detect the subtle signal that enabled us to implement the recharge mechanism proposed associated to the double-chamber model in agreement with recent petrological results.

The model presented in this study explains the uplift episode observed during the eruption of this volcano as a result of a magma flow from a deeper to a shallower reservoir. On this assumption, this model

could be used for forecasting the evolution of future eruptions, as it provides an estimate of the amount of magma that could be rising and how long it would take to reach a steady state. Our model can also offer a simple explanation for the differences in far-field and near-field displacement observed at many other volcanoes that not end in eruption and provides clues on how the recharge mechanism occurs in the different stages of the eruptive cycle including pre-eruptive, co-eruptive and post-eruptive stages.

## CRediT authorship contribution statement

C. Novoa Lizama: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. D. Remy: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. J.C. Baez: Writing – review & editing, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. A. Oyarzun: Writing – review & editing, Methodology, Investigation, Conceptualization. S. Bonvalot: Writing – review & editing, Resources, Funding acquisition, Conceptualization. A. Hooper: Writing – review & editing, Visualization, Resources, Investigation, Funding acquisition, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

We thank the European Space Agency and CNES for providing us Sentinel-1 data. ALOS-2 data from the ALOS satellite mission operated by the Japanese Aerospace Exploration Agency (JAXA) were used under the terms and conditions of the fourth ALOS 2 Research Announcement (project 1142). This work contributes to the International Research Network ‘Andes-FrenSZ’ (IRD, CNRS). We thank Valerie Cayol for guidance in using DEFVOLC and we are grateful to the Mesocentre Clermont Auvergne University (https://mesocentre.uca.fr/) for providing help and/or computing and/or storage resources. CNL and AH acknowledge support from the European Research Council (ERC) through the EU Horizon 2020 project DEEPVOLC (Grant 866085) and the UK Natural Environment Research Council (NERC) through the Centre for the Observation and modeling of Earthquakes, Volcanoes and Tectonics (COMET, http://comet.nerc.ac.uk). JCB thank the support of the project FONDECYT Regular N°1240501, ANID, Chile. Finally, we thank the editor Sonia Calvari and the two anonymous reviewers for their helpful and constructive suggestions that greatly improved this study.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jvolgeores.2024.108253.

## Data availability

No data was used for the research described in the article.

## References

Astor, A., Boixart, G., Folguera, A., Battaglia, M., 2022. Volcanic unrest at Nevados de Chillán (Southern Andean Volcanic Zone) from January 2019 to November 2020, imaged by InSAR. J. Volcanol. Geotherm. Res. 427, 107568. https://doi.org/10.1016/j.jvolgeores.2022.107568.
Bato, M.G., Pinel, V., Yan, Y., Jouanne, F., Vandemeulebrouck, J., 2018. Possible deep connection between volcanic systems evidenced by sequential assimilation of geodetic data. Sci. Rep. 8, 11702. https://doi.org/10.1038/s41598-018-29811-x.

Benet, D., Costa, F., Pedreros, G., Cardona, C., 2021. The volcanic ash record of shallow magma intrusion and dome emplacement at Nevados de Chillán Volcanic complex, Chile. J. Volcanol. Geotherm. Res. 417, 107308. https://doi.org/10.1016/j.jvolgeores.2021.107308.

Biggs, J., Pritchard, M.E., 2017. Global volcano monitoring: what does it mean when volcanoes deform? Elements 13, 17-22. https://doi.org/10.2113/gselements.13.117.

Biggs, J., Ebmeier, S.K., Aspinall, W.P., Lu, Z., Pritchard, M.E., Sparks, R.S.J., Mather, T. A., 2014. Global link between deformation and volcanic eruption quantified by satellite imagery. Nat. Commun. 5, 3471. https://doi.org/10.1038/ncomms4471.

Cardona, C., Gil-Cruz, F., Franco-Marín, L., San Martín, J., Valderrama, O., Lazo, J., Cartes, C., Morales, S., Hernández, E., Quijada, J., et al., 2021. Volcanic activity accompanying the emplacement of dacitic lava domes and effusion of lava flows at Nevados de Chillán volcanic complex Chilean Andes (2012 to 2020). J. Volcanol. Geotherm. Res. 736 (420), 107409.

Cayol, V., Cornet, F.H., 1998. Effects of topography on the interpretation of the deformation field of prominent volcanoes - application to Etna. Geophys. Res. Lett. 25 (1979-1982), 740. https://doi.org/10.1029/98GL51512.

Cembrano, J., Lara, L., 2009. The link between deformation and tectonics in the Southern Volcanic Zone of the Chilean Andes: a review. Tectonophysics 743 (471), 96-113.

Champenois, J., Pinel, V., Baize, S., Audin, L., Jomard, H., Hooper, A., Alvarado, A., Yepes, H., 2014. Large-scale inflation of Tungurahua volcano (Ecuador) revealed by persistent scatterers SAR interferometry. Geophys. Res. Lett. 41, 5821-5828. https://doi.org/10.1002/2014GL060956.

Christopher, T.E., Blundy, J., Cashman, K., Cole, P., Edmonds, M., Smith, P.J., Sparks, R. S.J., Stinton, A., 2015. Crustal-scale degassing due to magma system destabilization and magma-gas decoupling at Soufrière Hills Volcano, Montserrat. Geochem. Geophys. Geosyst. 16 (2797-757), 2811. https://doi.org/10.1002/2015GC005791.

Coppola, D., Laiolo, M., Lara, L.E., Cigolini, C., Orozco, G., 2016. The 2008 "silent" eruption of Nevados de Chillán (Chile) detected from space: effusive rates and trends from the MIROVA system. J. Volcanol. Geotherm. Res. 327, 322-329. https://doi.org/10.1016/j.jvolgeores.2016.08.016.

Del Negro, C., Currenti, G., Scandura, D., 2009. Temperature-dependent viscoelastic modeling of ground deformation: application to Etna volcano during the 1993-1997 inflation period. Phys. Earth Planet. Inter. 172, 299-309.

Dixon, H.J., Murphy, M.D., Sparks, S.J., Chávez, R., Naranjo, J.A., Dunkley, P.N., Young, S.R., Gilbert, J.S., Pringle, M.R., 1999. The geology of Nevados de Chillán Volcano, Chile. Rev. Geol. Chile 26, 227-253.

Dragoni, M., Magnanensi, C., 1989. Displacement and stress produced by a pressurized, spherical magma chamber, surrounded by a viscoelastic shell. Phys. Earth Planet. Inter. 56, 316-328. https://doi.org/10.1016/0031-9201(89)90166-0.

Dunkley, P., Young, S., 2000. Volcanic Hazard Mapping for Development Planning.

Farías, C., Lupi, M., Fuchs, F., Miller, S.A., 2014. Seismic activity of the Nevados de Chillán volcanic complex after the 2010 mw8. 8 Maule, Chile. Earthquake. J. Volcanol. Geotherm. Res. 283, 116-126.

Global Volcanism Program, 2023. Report on Nevados de Chillán (Chile). In: Sennert, S. (Ed.), Weekly Volcanic Activity Report, 4 January-10 January 781 2023.

González-Vidal, D., Obermann, A., Tassara, A., Bataille, K., Lupi, M., 2018. Crustal model of the southern central andes derived from ambient seismic noise rayleigh-wave tomography. Tectonophysics 744, 215-226.

Got, J.L., Peltier, A., Staudacher, T., Kowalski, P., Boissier, P., 2013. Edifice strength and magma transfer modulation at Piton de la Fournaise volcano. J. Geophys. Res. Solid Earth 118 (5040-5057), 791. https://doi.org/10.1002/jgrb.50350.

Heap, M.J., Villeneuve, M., Albino, F., Farquharson, J.I., Brothelande, E., Amelung, F., Baud, P., 2020. Towards more realistic values of elastic moduli for volcano modelling. J. Volcanol. Geotherm. Res. 390, 791.

Kositsky, A.P., Avouac, J.P., 2010. Inverting geodetic time series with a principal component analysis-based inversion method. J. Geophys. Res. Solid Earth 115. https://doi.org/10.1029/2009JB006535.

López-Escobar, L., Cembrano, J., Moreno, H., 1995. Geochemistry and tectonics of the Chilean Southern Andes basaltic Quaternary volcanism (37-46 S). Andean Geol. 22, 219-234.

Lupi, M., Trippanera, D., Gonzalez, D., D'amico, S., Acocella, V., Cabello, C., Stef, M.M., Tassara, A., 2020. Transient tectonic regimes imposed by megathrust earthquakes and the growth of NW-trending volcanic systems in the Southern Andes. Tectonophysics 774 (228204), 803. https://doi.org/10.1016/j.tecto.2019.228204.

Mattioli, G.S., Dixon, T.H., Farina, F., Howell, E.S., Jansma, P.E., Smith, A.L., 1998. GPS measurement of surface deformation around Soufrière Hills Volcano, Montserrat from October 1995 to July 1996. Geophys. Res. Lett. 25, 3417-3420. https://doi.org/10.1029/98GL00931.

Mee, K., Gilbert, J.S., McGarvie, D.W., Naranjo, J.A., Pringle, M.S., 2009. Palaeoenvironment reconstruction, volcanic evolution and geochronology of the Cerro Blanco subcomplex, Nevados de Chillán volcanic complex, central Chile. Bull. Volcanol. 71, 933-952. https://doi.org/10.1007/s00445815009-0277-7.

Moussallam, Y., Bani, P., Schipper, C.I., Cardona, C., Franco, L., Barnie, T., Amigo, A., Curtis, A., Peters, N., Aiuppa, A., et al., 2018. Unrest at 818 the Nevados de Chillán volcanic complex: a failed or yet to unfold magmatic 819 eruption? Volcanica 1, 19-32.

Moussallam, Y., Barnie, T., Amigo, A., Kelfoun, K., Flores, F., Franco, L., Cardona, C., Cordova, L., Toloza, V., 2021. Monitoring and forecasting hazards from a slow growing lava dome using aerial imagery, tri-stereo Pleiades-1A/B imagery and PDC numerical simulation. Earth Planet. Sci. Lett. 564, 116906. https://doi.org/10.1016/j.epsl.2021.116906.

Naranjo, Jose, A., Lara, L.E., 2004. August-September 2003 small vulcanian eruption at the Nevados de Chillán Volcanic Complex (36°50'S), Southern Andes (Chile). Rev. Geol. Chile 31, 359-366.

Newman, A.V., Dixon, T.H., Gourmelen, N., 2006. A four-dimensional viscoelastic deformation model for Long Valley Caldera, California, between 1995 and 2000. J. Volcanol. Geotherm. Res. 150 (832), 244-269.

Novoa, C., Remy, D., Gerbault, M., Baez, J.C., Tassara, A., Cordova, L., Cardona, C., Granger, M., Bonvalot, S., Delgado, F., 2019. Viscoelastic relaxation: a mechanism to explain the decennial large surface displacements at the Laguna del Maule silicic volcanic complex. Earth Planet. Sci. Lett. 521, 46-59. https://doi.org/10.1016/j.epsl.2019.06.005.

Novoa, C., Gerbault, M., Remy, D., Cembrano, J., Lara, L.E., Ruz-Ginouves, J., Tassara, A., Baez, J.C., Hassani, R., Bonvalot, S., Contreras-Arratia, R., 2022. The 2011 Cordón Cauille eruption triggered by slip on the Liquine-Ofqui fault system. Earth Planet. Sci. Lett. 583, 117386. https://doi.org/10.1016/j.epsl.2022.117386.

Orozco, G., Jara, G., Bertin, D., 2016. Peligros del Complejo Volcánico Nevados de Chillán: Región del Biobío. Servicio Nacional de Geología y Minería, Subdirección Nacional de Geología.

Oyarzun, A., Lara, L.E., Tassara, A., 2022. Decoding the plumbing system of Nevados de Chillán Volcanic complex, Southern Andes. J. Volcanol. Geotherm. Res. 422, 107455. https://doi.org/10.1016/j.jvolgeores.2021.107455.

Pedreros-Delgado, G., Contreras Vargas, M., San Martin, J., Cordova Varas, L., Bucarey, C., 2024. Análisis de actividad superficial sin precursors evidentes: El caso de la erupción (2016-2022) de Nevados de Chillán, chile. Cities on Volcanoes 12, COVID12 february 11-17.

Pritchard, M.E., Mather, T.A., McNutt, S.R., Delgado, F.J., Reath, K., 2019. Thoughts on the criteria to determine the origin of volcanic unrest as magmatic or non-magmatic. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 377, 20180008. https://doi.org/10.1098/rsta.2018.0008.

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., Bascou, P., 2014. A two-magma chamber model as a source of deformation at Grífsvötn Volcano, Iceland. J. Geophys. Res. Solid Earth 119, 4666-4683. https://doi.org/10.1002/2013JB010569.

Rodríguez-Molina, S., González, P.J., Charco, M., Negredo, A.M., Schmidt, D.A., 2021. Time-scales of inter eruptive volcano uplift signals: three sisters volcanic center, oregon (united states). Front. Earth Sci. 8, 577588.

Segall, P., 2016. Repressurization following eruption from a magma chamber with a viscoelastic aureole. J. Geophys. Res. Solid Earth 121, 8501-8522. https://doi.org/10.1002/2016JB013597.

Takeuchi, S., 2011. Preeruptive magma viscosity: an important measure of magma eruptibility. J. Geophys. Res. Solid Earth 116. https://doi.org/10.1029/2011JB008243.

Wasser, V.K., Lopez, T.M., Anderson, K.R., Izbekov, P.E., Freymueller, J.T., 2021. Multidisciplinary constraints on magma compressibility, the pre-eruptive exsolved volatile fraction, and the h2o/co2 molar ratio for the 2006 augustine eruption, alaska. Geochem. Geophys. Geosyst. 22, e2021GC009911.

Wicks, C., de la Llera, J.C., Lara, L.E., Lowenstein, J., 2011. The role of dyking and fault control in the rapid onset of eruption at Chaitén volcano, Chile. Nature 478, 374-377. https://doi.org/10.1038/nature10541.

Yip, S.T.H., Biggs, J., Edmonds, M., Liggins, P., Shorttle, O., 2022. Contrasting volcanic deformation in arc and ocean island settings due To 887 exsolution of magmatic water. Geochem. Geophys. Geosyst. 888 (23), e2022GC010387. https://doi.org/10.1029/2022GC010387.