Benefits and Limitations of the Growth Inversion Approach in Volcano Gravimetry Demonstrated on the Revisited 2004–2005 Tenerife Unrest

We review the current geoscientific knowledge of the volcanic unrest of 2004–2005 on Tenerife (Canary Islands) and revisit its gravimetric imprint. We revise the interpretation of the observed spatiotemporal (time-lapse) gravity changes accompanying the unrest by applying the Growth inversion approach based on model exploration and free geometry growing source bodies. We interpret the Growth solution, our new gravimetric model of the unrest, in the context of structural controls and the existing volcanological and geological knowledge of the central volcanic complex (CVC) of the island. Structural controls are inferred from the updated structural subsurface CVC density model obtained by our new Growth inversion of the available complete Bouguer anomalies (CBA data). Our gravimetric picture sees the unrest as a failed eruption, due to a stalled magma intrusion in the central position below the Teide–Pico Viejo stratocones, followed by upward and lateral migration of volcanic fluids reaching the aquifer and the SW end of the caldera wall. We thus classify the volcanic unrest of 2004–2005 as hybrid, in agreement with previous studies. The Growth inversion indicates that magma propagated along the boundary between the basaltic core of the island, the Boca Tauce volcanic body and the more permeable (less compacted) volcanic rocks with lower density. This gravimetric picture of the unrest provides new insights into the potential future reactivation of the volcanic system. Current geoscientific knowledge of the Tenerife volcanic unrest of 2004–2005 is reviewed New insights into the unrest are yielded by Growth inversion of observed time-lapse gravity changes Role of the freely adjustable inversion parameters in the Growth methodology is demonstrated Pros and cons of the Growth inversion approach in volcano gravimetric applications are illustrated Current geoscientific knowledge of the Tenerife volcanic unrest of 2004–2005 is reviewed New insights into the unrest are yielded by Growth inversion of observed time-lapse gravity changes Role of the freely adjustable inversion parameters in the Growth methodology is demonstrated Pros and cons of the Growth inversion approach in volcano gravimetric applications are illustrated


Introduction
Time-lapse gravity changes observed in volcanic areas during unrest can detect and characterize subsurface processes during volcanic unrest before the magma forces its way to the surface and triggers an eruption. Alternatively these data may help distinguish and characterize changes within hydrothermal systems. Spatiotemporal gravity changes are an observable signal produced by subsurface mass and/or density changes such as the migration of mobile masses (intrusion, injection, extrusion or withdrawal of magma or fluids), or changes in density of static masses, of magma residing within the plumbing system (due to physical or compositional changes induced by cooling or heating, fractional crystallization, mixing, degassing or vesiculation), as well as of the rock medium (due to compression, thermal expansion and displacement of density interfaces). Repeat (time-lapse) microgravity observations on benchmarks of a gravimetric network (Battaglia et al. 2008) may be complemented by continuous gravity recording at one or several stations to capture the faster intra-campaign changes (Williams-Jones et al. 2008;Carbone et al. 2020).
Spatiotemporal gravity changes are coupled with surface and inner deformations (e.g., Battaglia et al. 2008;Vajda et al. 2019), so joint analysis of both the gravity changes and the deformations can reduce any ambiguity in their interpretation. This can be done by joint analytical or numerical modelling (Battaglia et al. 2008). Observations of geodetic deformation are only sensitive to changes in the source pressure or volume (Rundle 1982;Fernández and Rundle 1994;Fernández et al. 1997Fernández et al. , 2006Fernández et al. 2021), but are blind to changes in source density and cannot therefore yield definite clues to the nature of the source process. This ambiguity can be resolved by observing changes in gravity which can constrain the mass of the intrusion (Rundle 1982;Fernández and Rundle 1994;Fernández et al. 2001a, b). The combination of deformation and 4D microgravity measurements can indicate the source density and better define the source process of the unrest (Fernández et al. 2001a, b;Battaglia et al. 2006;Camacho et al. 2011a). This is true if both the volume change generating the surface deformation and the mass change generating the gravity changes occur at the same location and in the same volumetric domain, when the joint modelling of deformations and gravity changes is of great benefit.
There are situations when significant gravity changes are observed accompanied by no measurable vertical surface displacements (e.g., Rymer et al. 1993;Fernández et al. 2001b;Carbone et al. 2003;Gottsmann et al. 2006). In other cases, the source of the surface deformation produces no significant gravity changes while the source of the gravity changes produces no significant deformation. The two sources thus differ in their location and must be sought as separate sources (e.g., Miller et al. 2017a;Vajda et al. 2021). In this situation, the spatiotemporal gravity changes properly corrected for the gravitational effect of the surface deformation (decoupled from the surface deformation) can be inverse-modelled as a stand-alone quantity (e.g., Vajda et al. 2021). When the two sources differ, the "density of the source", which points to the nature of the process, cannot be determined by combining the volume change obtained from deformation modelling and the mass change obtained from gravity change modelling; in this case, the nature of the process can be inferred from the location and spatial character of the mass change source in the context of the geological structure determined by independent studies (e.g., Miller et al. 2017b).
The gravitational signal of hydrological effects (changes in ground water level and soil moisture) is typically present in the gravity changes and masks the signal from volcanological sources. It is a major challenge to correct for the hydrological signal, as it would 1 3 strictly speaking require detailed 4D modelling involving near-surface 3D knowledge of the structure and the change in the hydrological situation over the time span of the gravity change data (e.g., Hemmings et al. 2016).
Here we focus on inverting and interpreting spatiotemporal gravity changes as a standalone quantity. In particular, we concentrate on inversion approaches that make no a priori assumptions about the number or shapes of the source bodies sought. We revisit the gravity changes of the Tenerife (Canary Islands) 2004-2005 unrest (Fig. 1) acquired and processed by Gottsmann et al. (2006). We apply a free geometry, model exploration and subsurface cell population inversion approach called Growth (Camacho et al. 2021b). Our objective is to analyse the benefits and limitations of the Growth inversion by analysing the Growth gravimetric model of the unrest in the context of the current geophysical knowledge of the unrest. We also aim to bring new insights into the unrest based on the Growth approach.

Review of geophysical interpretations of the 2004-2005 Tenerife unrest
In the spring of 2004, almost a century after the latest eruption on the volcanic island of Tenerife in 1909, a significant increase in the number of onshore seismic events marked the reawakening of the central volcanic complex (CVC). The increase in onshore seismicity, including five felt earthquakes, coincided with both an increase in the diffuse emission of carbon dioxide in the zone known as the Santiago Rift (Pérez et al. 2005) and increased fumarolic activity at the summit of Teide volcano . Gottsmann et al. (2006) carried out time-lapse microgravity measurements between May 2004 and July 2005, at 4 repeat campaigns (May 2004, July 2004, April 2005and July 2005. They observed spatiotemporal gravity changes which indicated that the reactivation was accompanied by an addition to subsurface mass, although without any widespread surface deformation . Gottsmann et al. (2006)  addition of 1.1 · 10 11 kg (with lower and upper bounds of 8.4 · 10 9 kg and 2.0 · 10 11 kg, respectively) using a Gaussian Quadrature integration over the area of positive gravity change. They considered three scenarios to explain the observed gravity changes: magma intrusion at depth, migration of hydrothermal fluids and a hybrid of both. They ruled out shallow dyke emplacement. They noted that a dyke emplacement at greater depth (a few km below the surface) into the Santiago Rift, perhaps recharging an existing reservoir, or in the form of magma injection into a conjugated fault system, could not be unambiguously excluded. Finally, they attributed these gravity changes to the upward (and lateral) migration of hydrothermal fluids released from a deep reservoir, triggered likely by magma injection into the Santiago rift zone at greater depth (prior to May 2004). They present a schematic model for the unrest (Fig. 2a, sketched after Fig. 4), in which they show a N-S vertical section crossing the Santiago ridge close to the west of Pico Viejo and passing through the western part of the Las Cañadas caldera. In their conceptual model (cf. Fig. 2a) of events between May 2004 and July 2005, the following events took place: (1) the likely injection of magma during the peak of onshore seismic activity 2 weeks before installation of the gravimetric network in May 2004, (2) the release of fluids or perturbation of existing hydrothermal system causing migration of fluids from NW to SE (May-July 2004) and later (July 2004-April 2005 and further into July 2005) along a N-S striking zone, (3) the upward migration of fluids along the upper surface of the high-density low-permeability Boca Tauce magmatic body situated beneath the western caldera (Ablay and Kearey 2000), (4) fluid migration into an overlying aquifer, located at a depth greater than 900 m beneath the caldera floor and thought to feed the Teide summit fumaroles (Araña et al. 2000). They note also that the western caldera boundary fault may act as a pathway for fluids to shallower depth.

estimated a mass
Based on seismic data from three antennas deployed in the Las Cañadas caldera between May and July 2004, Almendros et al. (2007) proposed a model for the unrest whereby a deep magma injection triggered a release of volatiles into an overlying aquifer. In their interpretation the magmatic reactivation started in April 2004 on the NW flank of the Teide at depths below 10 km b.s.l., related to an alleged basaltic magma chamber. The subsequent magma degassing triggered the rise of volatiles, which later reached the aquifer of the caldera (see their Fig. 15) and initiated the low-amplitude tremor on May 18 which lasted for several weeks.
The inversion of the Tenerife gravity data (complete Bouguer anomalies) by Gottsmann et al. (2008) provides evidence of a structural control of fluid flow beneath the CVC. Their model identifies pronounced low-density anomalies, which could indicate weaker permeable zones and preferred pathways for fluid migration.
Based on enhanced data processing and a new seismic catalogue, and by separating the volcano tectonic (VT) events into distinct families, Cerdeña Domínguez et al. (2011) have discovered the existence of two well-defined seismogenic zones related to the unrest: one located NW of the Teide-PV complex with mean depth 7 km b.s.l. and the other located on the SW border of Las Cañadas caldera with mean depth 4 km b.s.l. (Fig. 2b). They proposed a multiple source model for the 2004/2005 unrest and suggested a complex succession of deep-and shallow-seated magmatic intrusions beneath the CVC between 2001 and 2005. They reported alternating activity between the SW and NW seismogenic zones, which supposedly indicated an internal link between them, and maintained that this refutes the hypothesis of a gradual migration from N to S reported by Almendros et al. (2007) and . Although more than 95% of the seismic energy was released by the NW seismogenic zone (NWSZ), Cerdeña Domínguez et al. (2011) proposed a four-stage treelike magma intrusion model to explain the 2004 seismo-volcanic unrest. In their model 1 3 (see Fig. 2 and our sketch in Fig. 2b), which covers the time span 2001-2004, magma ascended from a greater depth at a central location below the Teide and was injected sideways into the NW and SW seismogenic zones, controlled by the rifts. The VT swarms observed in their seismic records were assumed to be caused by diking. The authors associate the NWSZ with the NW (Santiago) rift and the SWSZ with the prolongation of the NE (Dorsal) rift. The alternating seismic activity between the two zones gives reason for their hypothesis of a physical link between the two zones. The assumed physical link and the association of the VT events with diking gives in turn reason for their hypothesis on  Gottsmann et al. (2006) shown as the S-N vertical section at easting roughly 335,000 m, b a sketch of the tree-like magma model of Cerdeña Domínguez et al. (2011) inferred from seismic data, shown in lateral view, c a sketch of the gravimetric model of Prutkin et al. (2014) shown in lateral view the common central magma feeder (common intrusion origin) roughly aligned with the Teide-PV cones. The authors also quote the relation between basaltic intrusions and reactivation of phonolitic eruptions pointed out by . Cerdeña Domínguez et al. (2011) also bring up an argument worth considering. They liken the 1705-1706 rift eruptions to their proposed tree-like alternating-activity central-feeder magma model. At the beginning of 1705 a series of eruptions occurred in the NE rift (Siete Fuentes, Fasnia and Arafo) and approximately a year later, in May 1706, a new volcano appeared in the NW rift (Garachico) separated by more than 25 km from the eruptive vents of the previous year. They hypothesize that despite the spatial separation of these eruptions, the short time lag suggests the occurrence of a single eruptive process, although they admit that petrological evidence to confirm this hypothesis is lacking.
The spatiotemporal gravity changes of the unrest were revisited by Prutkin et al. (2014), who interpreted (see Fig. 2c) the 2004/2005 non-eruptive unrest also as hybrid (rather than purely magmatic or purely hydrothermal), caused by stalled intrusion of magma (mass injection estimated at 15 e + 10 kg) to a depth of about 6 km b.s.l. (deep source "Dp") roughly at the centre of the NW seismogenic zone (NWSZ) identified by Cerdeña Domínguez et al. (2011), about 5 km to the NW of the Teide-PV summits. The intrusion triggered the upward release of volcanic fluids (shallow sources "ShSW", "ShNW1" and "ShNW2" of mass estimated at 2 · 10 10 kg) through a network of permeable pathways or zones of least resistance into shallow-seated hydrothermal systems, causing the perturbation of the aquifers. Although the shallow SW source is vertically aligned with the centroid of the SW seismogenic zone (SWSZ) of Cerdeña Domínguez et al. (2011), it is not clear whether these fluids could have migrated upwards vertically as a result of magma injection into the SWSZ in January 2004 advocated by Cerdeña Domínguez et al. (2011). More likely they migrated laterally from the location of the shallow NW sources (cf. Fig. 2d), which correlates better with the fluids migration aspects of the models by Gottsmann et al. (2006) and Almendros et al. (2007). The gravimetric model of Prutkin et al. (2014) has not detected any mass addition at the locus of the SWSZ of Cerdeña Domínguez et al. (2011). Neither has it detected any mass addition at the trunk of their tree-like magma model, i.e. at the central location below Teide (compare Fig. 2c with Fig. 2b). The deep source ("Dp") representing the magma intrusion/injection ( Fig. 2c) correlates better with that of Almendros et al. (2007), although it is located further NW (also understanding that the depth reach of the gravimetric model goes only to about 7 km b.s.l.).
There are basically two views of the temporal evolution of the unrest. The first view adopts a north to south gradual migration of the hybrid (magma injection followed by released fluids) unrest sources (Gottsmann et al. 2006;Almendros et al. 2007;). This view is supported by temporal evolution of the shallow onshore seismicity. The anomalous inland seismicity was initially located beneath the northern and north-western flanks of the Teide-PV complex along the western portion of the Icod valley, but subsequently migrated towards the interior of the island along the western sector of Las Cañadas caldera ). This view is supported also by temporal evolution of the gravity high of the observed gravity changes. A marked gravity high (of 40 μGal amplitude) was found in the NW part of the covered area. With time, this high appeared to migrate in a southward direction, reaching the western part of the Las Cañadas caldera between July 2004 and April 2005. The gravity increase along the Santiago del Teide rift, noted between the first two campaigns (May and July of 2004), had disappeared by April 2005. At the same time, gravity increased significantly along the northern slopes of Teide, adding to the impression of a spatiotemporal migration of the causative source (Gottsmann et al. 2006;). In this view a deep intrusion of magma into the Santiago del Teide rift and a subsequent migration of hydrothermal fluids is favoured to explain the seismic and gravity data (Gottsmann et al. 2006;Almendros et al. 2007).
The second view, advocated by Cerdeña Domínguez et al. (2011), favours the alteration of unrest activity between the NW and the SW seismogenic zones, separated by 10 km, starting in January 2004, i.e. 3 months earlier than claimed in the previous studies. The authors propose a competing model, which concurs with previous geophysical and geochemical studies of the 2004 crisis. This model proposes a single magma intrusion affecting the central part of the island with lateral dikes driven by the rifts to the NW and SW.
Among other observations, summarized by , testifying that the unrest was volcanic in origin and nature, are increased (above usual background level) diffuse degassing carbon dioxide levels at the Santiago rift zone, and at the Teide summit crater, increased sulphur dioxide levels at the Teide crater, a fumarole at the Teide summit observed on 20 October 2004 and a new fracture (few tens of metres long) at the Orotava valley in the Benijos area opened on 5 December 2004 with associated gas emission.
When interpreting geophysical and geodetic data related to unrest, it is important to correlate the inversion results and the features of the model with structural geological or tectonic knowledge of the volcano, and with up-to-date knowledge of its plumbing system and preceding eruptive activity (Fernández et al. 2017(Fernández et al. , 2021. This facilitates the interpretation of the data observed on the unrest in the volcanological and geological context.

Geological structure of Tenerife
Any analysis of the subsurface geodynamic processes associated with the awakening or unrest of a volcanic system must be placed in the context of the knowledge of the potential pathways for ascending magma or the existing plumbing system. It is therefore essential to know its structural geology and tectonics. Tenerife, an intraplate volcanic island, has experienced an intriguing interplay between basaltic and phonolitic eruptive phases and events. Its construction is based on the basaltic shield followed by the central complex consisting of the Cañadas edifice and the twin stratovolcanoes of Teide and Pico Viejo placed within the Las Cañadas caldera. The Cañadas edifice is a composite volcano characterized by abundant explosive eruptions of highly evolved phonolitic magmas. The mostly effusive eruptive activity of Teide and Pico Viejo has evolved from basaltic to phonolitic. Below we review the knowledge of the island's geological structure as inferred from geological and geophysical data. A comprehensive summary is given in (e.g., Sainz-Maza Aparicio et al. 2019).
The pattern of evolution of Tenerife and the Canary archipelago is that of oceanic volcanic islands (Carracedo and Troll 2013). On Tenerife, first a shield volcano (or three shields) developed, formed mainly by basalts, trachybasalts and trachytes. The central shield was surrounded or accompanied by the Teno and Anaga shield volcanoes. Erosion and massive landslides devastated the northern flank of the shield (or shields). The shield phase was followed by post-erosive rejuvenation volcanism, mainly at the centre of the island, which formed the Las Cañadas volcano, which developed in three successive phases separated by large-scale flank collapses (Ancochea et al. 1999). Although Phase 1 was predominantly effusive and basaltic, the eruptions were more differentiated (trachybasalts and phonolites) and more explosive in Phases 2 and 3, when Plinian episodes produced pyroclastic falls and pyroclastic flows. Martí et al. (1997) proposed three main basaltic-tophonolitic cycles of development for the Las Cañadas volcano, each initiated with mafic or intermediate eruptions which then evolved towards phonolitic products. This succession of events seems to point to the simultaneous existence and interaction of rift zones and the felsic Las Cañadas volcano. The former are probably responsible for the basaltic (fissural) eruptions and the successive flank collapses. In this context, the development of the Las Cañadas Caldera and the Teide-Pico Viejo stratocones could represent the pinnacle of the latest cycle (Carracedo and Troll 2013).
These studies led to the recognition that Tenerife is essentially a basaltic shield structure on which a central phonolitic complex had developed. The basaltic shield is currently manifested on the surface at the three corners of the island: Teno, Anaga and Roque del Conde, while the central part is characterized by a high-density core covered by the products from the infill of the Las Cañadas caldera and from Teide-Pico Viejo (Ablay and Kearey 2000;Araña et al. 2000;Gottsmann et al. 2008). The uppermost part of this dense body is the old buried volcano of Boca Tauce. Other high-density regions were identified under the Teno, Anaga and Roque del Conde massifs. Las Cañadas has a negative gravity anomaly and its three main depressions, corresponding to the three overlapping caldera episodes (Martí et al. 1994;Martí and Gudmundsson 2000), have been clearly defined by magnetotelluric and electrical methods (Pous et al. 2002;Coppo et al. 2008;Villasante-Marcos et al. 2014).
Available geophysical data, supplemented with stratigraphic information from water galleries and borehole data, show that the caldera fill corresponds to deposits from Las Cañadas and at least 540 m of lavas and pyroclasts from Teide and Pico Viejo (Iribarren 2014). Observations suggest that the caldera collapse was as much as 1200 m. The three main flank collapse valleys, Güimar, La Orotava and Icod, were also filled with products from Las Cañadas, the rifts and Teide and Pico Viejo, except for the Icod valley, which has no products from Las Cañadas (Coello 1973;Iribarren 2014). Gottsmann et al. (2008) performed the inversion of the gravity data, the complete Bouguer anomalies (CBA), using the Growth inversion tool of Camacho et al. (2000Camacho et al. ( , 2002 for both homogeneous and stratified background densities. Their inversion results revealed that the subsurface anomalous density distribution model provides evidence of a structural control of fluid flow beneath the CVC. Their model identifies pronounced low-density anomalies which could indicate weaker permeable zones and preferred pathways for fluid migration. In a follow-up work, the same CBA data were inverted by Camacho et al. (2011b) using an upgraded version of the Growth software (GROWTH 2.0). The improved subsurface structural density model for Tenerife is presented in Fig. 6. Camacho et al. (2011c) performed a similar inversion of CBA data on Tenerife using a Growth approach modified to obtain subsurface structural density models dominated by subhorizontal layers delineated by irregular discontinuity interfaces. The modification consists of giving more weight in the solution to the subhorizontal layers rather than to isolated anomalous source bodies. This is achieved by taking into account an additional minimization condition, the proximity to prescribed horizontal interfaces, which is incorporated by including an additional weighting that is inversely proportional to the distance from the interface in the covariance matrix of the model parameters. This modified Growth approach unifies the inversion results based on isolated source bodies, basins with light infill and subhorizontal layered structures. The resulting structural density model for the subsurface of Tenerife is presented in Figs. 9 and 10.
A similar gravimetric study based on genetic algorithm inversion and focusing on the NW, central and S sectors of the island was conducted by Sainz-Maza Aparicio et al. (2019). They suggest the shallow character of the rift system and lateral migration pathways for magma along the rift, with ascending paths from a deeper zone located towards the east of this NW rift zone. They maintain that their model of subsurface density anomalies indicates the existence of several main channels in the north and south of the island that allow the ascent of deep basaltic magma to shallower zones of Tenerife.
Sainz-Maza Aparicio et al. (2019) used 3D gravity inversion procedures based on a genetic algorithm to image the internal structure of the shallowest part of the basaltic plumbing system of Tenerife and to explain how it controls magma transport and the location of basaltic volcanism. They compared the internal density structures of two zones that have recently been the site of basaltic volcanism: the Santiago del Teide rift area (NW rift), in which most of the cones are aligned parallel to the same structure, and the southern Tenerife area with an older monogenetic volcanic field, in which volcanic cones form several clusters and do not follow a single alignment. Their results suggest that the density anomalies observed below the NW rift zone are shallow and aligned parallel to the rift, which is consistent with a shallow rift system and a lateral migration of magma along it when the magma comes from a deeper zone located towards the east of this zone. The density anomalies in the south are also shallow but distributed along more than one alignment and are connected to similar but much deeper subvertical density anomalies. These results allegedly support the existence of several main channels in the north and south of the island that allow the ascent of deep basaltic magma to the shallower zones of Tenerife, where the magma may be either caught by the rift systems, with a lateral distribution and occasional eruptions, or continue its ascent to the surface using a secondary system of shallow fractures that dominate the southern part of Tenerife.

Plumbing system
A key aspect in understanding the internal structure of a volcano is the size and nature of its magmatic system. Volcanic hazards are related to the size and depth of the magma sources. The products of past eruptions show evidence of two styles of magma ascent on Tenerife: one associated with active rift basaltic volcanism and another associated with the explosive volcanism of phonolitic magmas. This implies that magma has been tapped from different regions at different depths.
Alkaline basaltic ridge magmatism implies a deeper mafic reservoir and different ascent pathways. From the geophysical perspective, neither magma storage system is clearly defined. De Barros et al. (2012) presented a seismic tomographic model of Tenerife which identifies two anomalous regions: one at ca. 7-9 km b.s.l. in the northern part of the island (UTM 3115-3140 km of the S-N cross section running through the Teide summit) and a shallower one at 1-4 km b.s.l. (in the W-E cross section directly below the Teide summit). They interpret the shallow anomaly as possibly pertinent to the alleged phonolitic magma storage feeding the Teide-PV complex and the deeper anomaly as potentially related to the basaltic feeder system. They mention two possibilities in their interpretation: these anomalous seismic regions represent either fractures filled by volatiles associated with the gasfilled upper part of the magma storage system or a complex set of dykes and sills filled with magma. Their deeper northern anomalous structure correlates well with the NW seismogenic zone of the 2004-2005 unrest of Cerdeña Domínguez et al. (2011). Piña-Varas et al. (2018) analysed the detectability of the magma storage regions using magnetotellurics (MT). MT studies have so far not detected any conductivity anomalies that could be clearly associated with the chambers. The authors report that shallower phonolitic chambers smaller than 3 × 3 × 1 km 3 could remain undetected by the existing MT surveys, and new data must be acquired to confirm or rule out their existence. They also suggest that no mafic reservoir is present above a depth of 8 km, and propose a model for the internal structure of Tenerife (see Fig. 5) by merging and comparing all the information derived from geological, geophysical and hydrogeochemical studies along a NW-SE profile.
To summarize, most of the geological, petrological and geophysical studies concur in locating the mafic reservoir associated with fissure eruptions in the northern part of the island at depths of between 5 and 14 km b.s.l., and the phonolitic chambers between 1 and 4 km b.s.l. The differences between phonolitic eruptions suggest that Teide-Pico Viejo is currently in the initial phase of magmatic evolution, producing less evolved and less explosive magmas stored in small chambers . Geological and geophysical studies (Canales 1997;Watts et al. 1997;Ablay et al. 1998;Neumann et al. 1999;Ablay and Kearey 2000;Dañobeitia and Canales 2000) point to periodic accumulation of basaltic magmas in reservoirs located at three structural discontinuities: the base of the lithosphere at 30 km b.s.l., the Moho discontinuity at 14-16 km b.s.l. and the ocean basement boundary of the Teide volcano at 7-8 km b.s.l. (Martí and Gudmundsson 2000).

Recent eruptive history
Knowledge of recent eruptions is also essential for understanding the unrest. Curiously enough, in 1492 Cristopher Columbus witnessed a volcanic eruption on Tenerife (now ascribed to the Boca Cangrejo vent), when he passed the islands on his voyage to discover the New World. Historical volcanic activity (Carracedo et al. 2007;Carracedo and Troll 2013) on the island is associated with vents on the Santiago (NW) rift (Boca Cangrejo in 1492, Montañas Negras in 1706, Chahorra in 1798 and El Chinyero in 1909) and the Dorsal (NE) rift (Fasnia in 1704, Siete Fuentes and Arafo in 1705). The 1706 Montañas Negras eruption was the most significant in recent history and destroyed the town and principal port of Garachico, which was the island's trade hub at the time, after forty days spitting out lava and volcanic matter. After the catastrophe, Tenerife's trade moved to Santa Cruz. The latest eruption to date was the rift-related basaltic strombolian eruption of 1909 at Chinyero.
The main concern, however, relates to the central active Teide-Pico Viejo volcanic complex, which poses serious hazards to the properties and population of the island and which has erupted several times during the last 5000 years, including a subplinian phonolitic eruption (Montaña Blanca) about 2000 years ago ). The last summit eruption from Teide occurred about 850 AD and produced the Lavas Negras that cover much of the flanks of the volcano. Historical activity within the Las Cañadas caldera, associated with the Teide and Pico Viejo stratovolcanoes, occurred in 1798 from the Narices del Teide (Teide's Nostrils, also known as Chahorra) on the western flank of Pico Viejo. The eruption was predominantly strombolian.

Methods applied to interpret the Tenerife 2004-2005 gravity changes
As for the methods applied to interpret the gravity changes of the unrest, Gottsmann et al. (2006) used inverse modelling for an infinite N-S striking horizontal cylinder (as a first order approximation), running just west of the Pico Viejo, to explore the scenario of fluids migration. Using a global optimization iterative method with various initial values for depth and radius, they found the source to be located at roughly 2 km below the surface. Prutkin et al. (2014) opted for a multisource inversion based on the so called method of local corrections. In their approach, they separated the field of the observed gravity changes into a shallow field and a deeper field. This required the interpolation of the scattered point data onto a regular grid and a removal of regional trend defined as a harmonic function, followed by a triple (upward-downward-upward) harmonic continuation of data. Then the shallow and the deeper fields were inverted separately by means of sources represented by mass line segments. The method usually involves an additional step in which the line segments are transformed into potato-shaped compact homogenous source bodies. This step was skipped and the sources of the unrest remained represented by the line segments. For the interpretation of their line segments, see Sect. 2 and Fig. 2c.
Here we use the Growth inversion approach, which belongs to the class of inversion methods that make no a priori assumptions about the number and shape of the sources (Camacho et al. 1997(Camacho et al. , 2002(Camacho et al. , 2011a(Camacho et al. , 2011b(Camacho et al. , 2020(Camacho et al. , 2021aCannavó et al. 2015;Fernández et al. 2021). It seeks the sources as being represented by a free aggregation of subsurface rectangular cells (prisms) filled in by a userprescribed constant differential density value or alternatively by a set of several density values (with a prescribed step) that may vary from cell to cell. The inversion is based on an automated iterative model exploration process using growing source bodies due to stepwise cell population.
The Growth inversion approach was originally developed (Camacho et al. 1997) for inverting complete Bouguer gravity anomalies (CBA) given on the topographic surface. The topographically corrected (Bouguer) gravity anomalies/disturbances represent the gravitational effect of subsurface anomalous density distribution (Vajda et al. 2006(Vajda et al. , 2020. The inversion approach underwent several improvements and refinements over time leading to the GROWTH-2 and GROWTH-3 versions (Camacho et al. , 2002(Camacho et al. , 2011b(Camacho et al. , 2021a and was applied in structural studies in several regions (e.g., ).
Here we use the most recent version of GROWTH to invert the Tenerife CBA data to produce a structural density model of the island to provide a background information on structural controls in the context of which we interpret the subsurface density changes obtained by Growth inverting the spatiotemporal gravity changes of the unrest.
The same approach turned out to be applicable also to the inversion of time-lapse gravity changes (Camacho et al. 2021b) or to the joint inversion of gravity changes and elevation changes (surface deformation) at benchmarks (Camacho et al. 2011a generated by mass changes and pressure changes, respectively. A similar approach was developed for inverting surface deformation fields (such as those obtained from satellite DIn-SAR data), based on seeking not only pressure sources but also free geometry dislocation sources, namely dip slip, strike slip and tensile sources (Camacho et al. 2020;Fernández et al. 2021). In the case of inverting spatiotemporal (time-lapse) gravity changes using the GROWTH-dg software application (Camacho et al. 2021b), the sources represent subsurface spatiotemporal (time-lapse) density changes that will be called simply "differential density" in the sequel, understanding they represent a subsurface density change distribution over a time interval.
The application of the Growth approach to spatiotemporal gravity changes raises new challenges. When inverting observed gravity changes, one commonly deals with undersampling of the information on the sources due to a limited number of observation points, which is contrary to inverting CBA data that are typically much more numerous. Also, while inverting the CBA data implies working with data at the level of tens of mGals, inverting gravity changes means dealing with data typically at the level of a few tens of μGals with observational error at the level of around 15 μGal, i.e. data with a much poorer signal-to-noise ratio. Due to the limited amount of input gravity data, the volcano gravimetric studies usually opt for inversions based on assumed simple-geometry source bodies in order to tackle an over-determined inverse problem, where the number of source parameters sought is smaller than the number of input gravity data. This approach is more suitable if there is external evidence or a justified reason for choosing the number of sources to seek (typically one to three) and their type (shapes/geometries), such as spheres and ellipsoids for magma chambers, prisms for sills and dikes and cylinders for chimney feeders.
On the other hand, "free geometry" inversion methods that do not rely on pre-constraining the number or shape of the sources have the advantage of not discriminating the sources beforehand. They can identify the subsurface locations of the most significant density changes, indicate the nature of the sources and shed light on the subsurface geodynamic process. We compared the two approaches in a case study at the Laguna del Maule volcanic field, Chile (Vajda et al. 2021), and found that the Growth approach may offer additional insights compared to the simple source analytical optimization inversion methods.

Input gravity data
The input gravity data for the Growth (or any gravimetric) inversion of gravity changes must correspond to the gravitational effect of the sources of temporal mass and/or density changes sought. By mass changes, we mean the addition or removal of mobile masses to/ from the subsurface volumetric domain of the model. Temporal density changes refer to the changes in density of the masses originally present within the subsurface model domain.
The observed gravity changes must therefore first be corrected for all signal components other than those associated with the geodynamic process under study, such as atmospheric and tidal effects, instrumental drift, hydrological effects (e.g., Jousset et al. 2000;Battaglia et al. 2008;Hemmings et al. 2016;Carbone et al. 2017;Fernández et al. 2017;Van Camp et al. 2017) and surface mass changes (e.g., Jousset et al. 2000). If surface deformation (uplift or subsidence) accompanies the observed gravity changes, they must be properly corrected for the gravitational effect of the surface deformation. This correction is given in terms of the deformation-induced topographic effect (Vajda et al. 2019) and is evaluated by Newtonian volumetric integration, unless it can be accurately estimated by a planar Bouguer or normal free air approximation (e.g., Berrino et al. 1984Berrino et al. , 1992Rymer 1994;Battaglia et al. 2008). The gravity changes are also affected by the gravitational effect of subsurface deformations (e.g., Rundle 1982;Fernández and Rundle 1994;Bonafede and Mazzanti 1998;Jousset and Okada 1999;Battaglia et al. 2006Battaglia et al. , 2008Vajda et al. 2019). However, unlike the surface vertical displacements, these deformations are not directly observable, so their effect cannot be computed and applied as a correction to the gravity changes, but only estimated or modelled. The application of all the above corrections to the gravity changes results in the residual gravity changes that are used as input data in the Growth inversion process.
We apply the inversion approach using the GROWTH-dg tool to the 4D microgravity data observed in the CVC of Tenerife during the seismo-volcanic unrest of 2004-2005, observed by Gottsmann et al. (2006). Specifically we focus on cumulative gravity changes observed between May 2004 (2 weeks after the start of the anomalous onshore seismic activity of the unrest) and July 2005. No significant (above 3 cm accuracy threshold) widespread surface deformation accompanying the gravity changes was observed (Gottsmann et al. 2006;Fernández et al. 2009). Consequently, no deformation-induced topographic correction (free air or Bouguer-corrected free air) was applied to the observed gravity changes.
A bulk gravity increase was recorded (Gottsmann et al. 2006) between May 2004 and July 2005 across 14 gravity benchmarks in the joint ground deformation/microgravity network (Fig. 1). No correction for the gravitational effect of surface deformation (e.g., Berrino et al. 1984Berrino et al. , 1992Battaglia et al. 2008;Carbone et al. 2017;Fernández et al. 2017;Van Camp et al. 2017;Vajda et al. 2019) was applied due to the lack of elevation changes observed at benchmarks above the threshold of 3 cm . Hydrological correction (Battaglia et al. 2003;Hemmings et al. 2016) was applied at several benchmarks based on data about changes in the water table level (Gottsmann et al. 2006). Tidal corrections, instrumental drift correction and adjustment of redundant observations were applied. The precision of the resulting gravity changes was estimated (Gottsmann et al. 2006) at 15 μGal (1 μGal = 10 -8 m/s 2 ).

Multitude of Growth solutions
We revisit the gravimetric signature of the Tenerife 2004/2005 unrest with a twofold aim. First, by applying the new GROWTH-dg tool to invert the spatiotemporal gravity changes of the unrest, we hope to contribute new additional insights to the nature of the unrest and the possible pathways for magma ascent in the CVC of Tenerife. Second, we aim to demonstrate the applicability, functionality and benefits of the Growth inversion approach in volcano gravimetric studies. We therefore set out to show the advantages and limitations of the Growth inversion methodology when applied to sparse time-lapse gravity changes with a low signal-to-noise ratio (in our case 40 vs 15 μGal).
We run the Growth inversion process on the point time-lapse gravity changes, given at benchmarks on the topographic surface, described above and shown (as interpolated) in Fig. 1. The GROWTH-dg tool has several inversion parameters that can be freely adjusted by the user (within certain limits) in each inversion run. The selection of the values for these parameters has an impact on the resulting inversion model and on the fit to the input gravity data. Since the input gravity data-the spatiotemporal gravity changes-are sparse and noisy, there will be a suite of diverse inversion models that produce a more or less acceptable fit to the data. These models cannot be discriminated based purely on the gravity data. Although this may appear as a drawback, it is in fact an advantage. It is preferable to have a multitude of models to choose from based on external knowledge, rather than having one particular model and wonder how realistic it might be. The admissible models can be discriminated based on external cognition, such as structural geological or tectonic information, volcanological considerations, independent geophysical data and their interpretation, eruptive products indicating magma pathways and eruptive styles, seismic unrest swarms or tremor location.
The user of the GROWTH-dg tool must be aware of the fact that the sparsity and poor signal-to-noise ratio of the input gravity data allows a relatively large suite of diverse inversion models to be obtained that almost equally well satisfy the input data, considering their error bars. Running the inversion with several values of the inversion parameters helps the user develop a feeling for the nature and spatial characteristics of the geodynamic process responsible for the 4D microgravity data observed. Below we show these inherent properties of the Growth inversion methodology by running the inversion on the Tenerife gravity changes with several different values of the individual inversion parameters, which produces a suite of admissible (acceptable based on the fit to the gravity data) yet diverse inversion models. We then select the model which in our view best fits both the gravity data and the available constraining information.
We demonstrate the sensitivity of the Growth inversion solutions (models), obtained for the 2004-2005 Tenerife unrest gravity changes, to the values of the user-preset values of the individual tunable inversion parameters in the online Supplement (Sect. 2). Here we show and discuss only several of those inversion models.
The inversion with a relatively high balance factor (λ = 35) and differential density of 10 kg/m 3 produces a simple compact model (Fig. 3) with one source body centred under Fig. 3 Inversion model obtained with λ = 35 and differential density of 10 kg/m 3 the NW slope of the twin stratovolcanoes at a depth of 6.7 km below sea level (b.s.l.), rising to a depth of about 5 km b.s.l., apparently from greater depths (limited by the maximum depth reach of the model domain at 8 km b.s.l.). The mass addition of the source body is 1.7 · 10 11 kg. This model has a somewhat poorer data fit (with misfit r.m.s. of 18 μGal and max 37 μGal). This solution is considered oversimplified due to the relatively high balance factor. However, such a solution is useful, because it indicates at which location the most significant mass/density changes may take place. This source body could be interpreted as a volumetric domain affected by magma intrusion/injection (swarm of thin dykes or dykes and sills, or eventually an injected conjugated faults system) causing the bulk density increase of 10 kg/m 3 on average over the entire source body volume.
The inversion solution with a smaller balance factor (λ = 20) and the same differential density (10 kg/m 3 ) is shown and discussed in the Supplement (Fig. S4). This model now consists of three source bodies (see Fig. S4 in the Supplement): the central, the southern (SW) and the northern (NW) bodies. For reasons discussed in detail in the Supplement, the northern source has to be considered as fictitious, an artefact of the inversion.
The input gravity data in volcano gravimetry (the time-lapse gravity changes) are often very sparse. For that reason, some of the data points may represent only local uncorrelated gravity signal, in other words, may bear information on site effects only. To suppress the contribution of such data points into the adjustment which aims to target sources other than those causing site effects, Growth offers the iterative reweighting functionality (inversion parameter B). We present two inversion solutions obtained with two different intensities of applied iterative reweighting. Figure 4 presents the model for B = 1.8 (with misfit r.m.s. of 11 μGal and mass injection 3.0 · 10 11 kg). Figure 5 shows the model for yet a stronger reweighting with B = 1.2 (with misfit r.m.s. of 8 μGal and mass injection 2.6 · 10 11 kg). The reweighting feature practically eliminates the contribution of two gravity benchmarks to the adjustment, which causes the model to change (compared to that of Fig. 3). Now, the source body is vertically elongated, reaching a depth level of 1 km b.s.l. (Fig. 4), or even more elongated with its fingers reaching a level of 1 km a.s.l. (Fig. 5). These two models produce a good fit, albeit possibly slightly too tight (referred to as "over-fitting"). The upper part of these two model source bodies reaches the depth level of the hypothetical, assumed based on geochemical/petrological analyses, phonolytic magma chamber/s. These two sources could also be interpreted as a volumetric domain affected by magma intrusion/injection, releasing volatiles in their upper parts, which could correspond with the gas emissions and fumarolic activity at the Teide summit.
In Fig. 6 we show a model in which the depth weighting is applied also, obtained with the following parameters: λ = 20, target density value 10 kg/m 3 (one level), B = 2.2 and D = 1. The model fit is good (with misfit r.m.s. of 11 μGal). The injected mass is estimated at 4.2 · 10 11 kg. This model seems to suggest that the northern source body is physically connected with the central body at the depth level of 7 km b.s.l. and deeper. However, such an impression appears illusory for several reasons. As explained above and in the Supplement, the northern source is to be considered as fictitious, as an inversion artefact only. The physical link appearing to connect the central and northern bodies is to be considered as fictitious also, produced by the depth weighting as an artefact. Higher values of the D parameter tend to produce extra flat bodies at the bottom boundary of the model domain, as demonstrated by synthetic case studies (work submitted for publication). Based on our experience from synthetic modelling, we do not recommend to use the depth weighting functionality when interpreting spatiotemporal gravity changes.

Interpretation based on preferred Growth model
To select the best candidate, out of the suite of Growth solutions, to shed light on the subsurface process, we consider not only the gravity data fit of the model, but also the geological and volcanological context. We use a structural density model of Tenerife obtained by a Growth inversion of the CBA data to spatially correlate the source bodies of mass changes, found by Growth inverse modelling the spatiotemporal gravity changes of the unrest, with the geological structure (see Figs. 7 and 8, cf. also Gottsmann et al. 2008).
All the Canary Islands include large high-density bodies composed of thousands of dykes, cooled chambers and risen basal complexes (Fernández et al. , 2021. In Tenerife this large core lies under the centre and south and suggests the volcanic sources that have formed the island, shown in shades of red in Figs. 7 and 8 (positive density contrast). It is a dense compact structure, perhaps still subjected to a slight isostatic readjustment, a subduction (Fernández et al. 2009) that is too hard to be penetrated by fluids and requiring possible rising magma to bypass it. Lower-density areas (shades of blue) correspond to unconsolidated deposits of volcanic and sedimentary materials. These areas are more likely to be penetrated by magma and/or volatiles.
The sparsely observed 4D microgravity data in volcanic areas often represent signals from uncorrelated shallow local/site effects, even if acquired and processed with due care. We illustrate this point on the gravity benchmark in the NW part located farthest west, which is the only one with a negative gravity change, and is thus isolated and not correlated with neighbouring gravity benchmarks. This benchmark coincides with a known local CBA minimum corresponding to low-density deposits, which likely houses a local aquifer. The drop in the water level in this aquifer may be the origin of the negative gravity change observed. The site also coincides with localized ground subsidence detected using SAR techniques interpreted as being due to the drop in the water level in the aquifer caused by water extraction (Fernández et al. 2005(Fernández et al. , 2009Gottsmann et al. 2008). This site effect may be present to some degree at other gravity benchmarks, see (Gottsmann et al. 2006) for a discussion of the hydrological effects relevant to the repeat gravimetric surveys covering the unrest.
All remaining gravity benchmarks exhibit a positive gravity change and contain a common correlated signal component which reaches its peak at the central benchmark. The resulting Growth inversion model (Figs. 7 and 8), obtained with a relatively high balance factor value (λ = 25) and showing only correlated components, indicates that the cause of this small positive gravity change is a vertically elongated positive deep source body (shown in purple in Figs. 7 and 8) just north (about 1.5 km) of the Teide-PV complex, near the boundary between the high-density (lower-permeability) island core body (Boca Tauce) and the low-density (weaker) rock medium, residing in the more permeable part.
Considering the location and shape of this source body, and taking into account the fact that the gravity changes were associated with seismic manifestation (Almendros et al. 2007;Cerdeña Domínguez et al. 2011) and temporal geochemical anomalies (gas emission changes), it is reasonable to interpret this body as a stalled magmatic intrusion that tapped a deeper reservoir beneath the centre of the island. Its ascent path appears to be tangential to the high-density core structure. Magma was intruded along a more permeable pathway due to the less compacted rock environment at the margin of the harder high-density Boca Tauce body, as shown in Figs. 7 and 8, and schematically also in Fig. 9. This intrusion released volatiles that migrated to shallower depths as hydrothermal fluids (at depths between sea level and 1 km a.s.l.) and disturbed the aquifers, which enabled lateral migration to the shallow source body at the SW end of the caldera wall as shown in Fig. 10.
Since our gravity changes span a time interval May 2004 to July 2005, we are not able to judge whether these migrated fluids (the SW source 1 km a.s.l.) were vertically

Benefits and limitations of the Growth inversion as applied to gravity changes
It is important for the user to realize that gravimetric inversion is a strongly ambiguous endeavour, especially when working with a small number of scattered and imprecise microgravity data. A whole spectrum of quite diverse inversion models can be obtained, all of them fitting (more or less properly) the input gravity data, as shown in Sect. 4.2 (and mostly in the Supplement) by our inversion models for the gravity changes of the Tenerife unrest. These models cannot be distinguished and interpreted based purely on the gravity data. External knowledge of the structural features and/or independent geophysical data and their interpretation (which are not themselves inherently free from ambiguities either) can be considered jointly to be able to select a proper inversion model from the multitude and interpret the resulting temporal subsurface density changes.
An additional challenge in inverting and interpreting spatiotemporal gravity changes using the Growth approach is that these microgravity data are very sparse and may largely represent a signal of uncorrelated local/site effects no matter how precisely and carefully the data were acquired in the field and processed.
The user should also understand that Growth cannot correctly detect the volume and differential density of sought sources. The differential density is preselected by the user to run the inversion process, not determined by the inversion. Moreover, the Growth approach naturally operates with small density contrasts that may be unrealistic with respect to the real sources generating the gravity data. This is partly due to the size of the cells of the subsurface partition. As a consequence the volumes of the found source bodies will be proportionally overestimated. This is particularly valid in case of thin elongated sources with high differential density such as dykes, sills, inclined sheets and chimneys for two reasons. One is that the average cell size of the subsurface partition is dictated by the spatial distribution of the gravity data points, which is typically at least a few hundred metres, while one or two dimensions of the source may be much smaller than the average cell size (way below the resolution capability). The other is that the Growth procedure must operate with small differential density values (preset by the user as an inversion parameter), in order to correctly identify the location, orientation and relative shape (geometry and hence the nature) of the sources. The Growth inversion approach senses the thin elongated sources as much thicker and much more voluminous, with a much smaller, bulk differential density, as demonstrated by synthetic case studies (work submitted for publication).
However, the mass of the injection/intrusion, as a product of the (found) source volume and its (preselected small) differential density, is correctly determined (in case of optimally chosen balance factor value, which means no over-fitting and no model oversimplification). This is an important aspect of the Growth inversion, since the mass of an intrusion/ injection is a key parameter in volcano gravimetry. In the case of magmatic intrusions, the source volume of small differential density may be interpreted as a swarm of thin dykes (or dykes and sills) or as filling a system of faults and fractures, whereby the differential density would represent a bulk change over the entire volume rather than the density of the injected fluid (magma). As a consequence, the Growth inversion cannot unequivocally distinguish between sources of injected magma and sources of brines filling porous medium.
In case of sources that are spatially correlated with existing or assumed magma chambers or reservoirs, the small density changes can be interpreted as being caused by physical or compositional changes in the stored or rejuvenated magma (cooling/heating, degassing/ vesiculation, mixing, fractional crystallization, etc.) In sources located in porous rock media the small density sources can be interpreted as migrated liquids and a change in the pore content.
The GROWTH-dg tool for microgravity data inversion offers numerous advantages thanks to its individual and tuneable inversion parameters. The user is advised to run numerous inversions with several values of each individual inversion parameter, sensitively varied during a trial-and-error process. The multitude of the inversion solutions obtained can then be analysed and validated against the independent structural geological, geophysical and volcanological information available for the given site and volcano.
Another continuing challenge when interpreting the subsurface temporal density change sources found by the Growth inversion is that these sources may represent both mass changes (mobile masses) and density changes (of masses already present), and cannot be easily distinguished without additional external information.

Lessons learned from revisiting the 2004 Tenerife unrest by the Growth inversion approach
Our preferred (best choice) Growth inversion model representing our new gravimetric picture of the 2004 unrest, and its interpretation, is presented in Sect. 5. In agreement with previous interpretations we conclude on a stalled intrusion (failed eruption) triggering migration of volcanic fluids. Thus we classify the volcanic unrest as hybrid. Since our inversion used the cumulative gravity changes observed between May 2004 and July 2005, our analysis cannot trace temporal evolution of the unrest over the given period, hence cannot contribute into the debate whether the gradual N to S migration perception of the unrest (Gottsmann et al. 2006;) or the alternating dyke propagation from the central intrusion into the NWSZ (within the Santiago rift zone) and SWSZ is preferable. Comparing our results to the seismic data-derived tree-like central magma intrusion model of Cerdeña Domínguez et al. (2011), we observe a correlation of our model with the proposed central magma intrusion (the trunk of the model). Our analysis detects the central magma feeding pathway as offset slightly (roughly 2 km) to the NW from the Teide summit, tracing the boundary between the higher-density (assumed lower-permeability) and lower-density (assumed higher permeability) rock environments. The structural density model of the CVC of Tenerife was obtained by us here using an upgraded Growth inversion of the CBA gravity data. Although the bottom reach of our model is 7 km b.s.l., there are clear signs that the central magma intrusion comes (subvertically) from greater depths. As for the upper reach of the central intrusion, we cannot unambiguously conclude based on our gravimetric inversion at what depth the magma intrusion stalled, from where its released fluids continued the upward ascent. What appears clear is the potential of this pathway to reach and affect the location of the assumed phonolitic magma chamber system (if it exists).
Our gravimetric inversion detected a mass addition source at the SW end of the caldera wall, at a depth span between sea level and 1 km a.s.l., interpreted by us as arrival of migrated hydrothermal liquids. This source correlates perfectly with the shallow SW source found by the gravimetric inversion of Prutkin et al. (2014) using a very different inversion technique and interpreted likewise. This shallow SW source lies straight above the locus (at 4 km b.s.l.) of the SWSZ inferred by Cerdeña Domínguez et al. (2011). However, neither of the two gravimetric methods detected a mass addition source at the locus of the SWSZ. Thus, if the fluids of the shallow SW source were to had been vertically released by magma intrusion into the SWSZ, it would have to had happed prior to May 2004, or the mass of the intrusion was too low to be gravimetrically detected at the depth of 4 km b.s.l. (at the threshold of 15 μGal on the surface). Another, more likely possibility seems to be that the fluids migrated from the central intrusion first vertically and at a shallower depth (via the aquifer) laterally to their final position.
As for the NWSZ and diking associated with it proposed by Cerdeña Domínguez et al. (2011), the locus of the NWSZ dwells at 7 km b.s.l., which is at the bottom boundary of our model and thus beyond the reach of the "trustworthy senses" of the gravimetric network. We did not detect gravimetrically any mass addition in the volumetric domain of the NWSZ above the level of 7 km b.s.l. Although the Growth models, obtained with balance factor values around 20 (or smaller), detect a source of mass addition in the north, in the Icod valley very close to the shore, with a locus at ca. 2.5 km b.s.l. and vertical span from sea level to roughly 5 km b.s.l. (cf. the online Supplement), this source is beyond the reach of the "trustworthy senses" of the gravimetric network, and is to be considered as fictitious (as inversion artefact) for reasons detailed in Sect. 4.2 and the Supplement. Moreover, this source lies well beyond the NWSZ, so it lacks support from seismic data.
Comparing our new gravimetric picture of the unrest with that of Prutkin et al. (2014), who used the same gravity data but a different inversion method, these two agree with seeing the unrest as hybrid and due to a stalled magma intrusion (failed eruption). There is an offset though in the horizontal position of the magma intruding from depth. The deep source (at roughly 6 km b.s.l.) of Prutkin et al. (2014) correlates with the locus of the NWSZ, while our central intrusion sources lies closer to the Teide-PV summits. There is agreement between our shallow SW fluids source and their shallow SW source interpreted likewise. However, the shallow NW sources found by Prutkin et al. (2014) are missing in our new gravimetric picture. Our gravimetric model of the unrest seems to correlate with the concept of the unrest proposed by Almendros et al. (2007).

Conclusions
The Growth inversion approach, with its un-preconstrained (free geometry) sources and model exploration with subsurface cells population and growing source bodies, and its software realization GROWTH-dg (Camacho et al. 2021a, b), has been shown to be useful in volcanogravimetric interpretation thanks to its versatility. It makes use of several tuneable inversion parameters to produce a multitude of more or less (based on the fit to the observed 4D microgravity data) admissible solutions (inversion models) that can be discriminated and interpreted when correlated with external knowledge, such as results from other independent geophysical methods and external geological and/or volcanological information.
We have revisited the gravimetric picture of the 2004-2005 Tenerife unrest by applying the Growth inversion approach to the spatiotemporal gravity changes accompanying the unrest. We interpret the unrest as hybrid and due to a stalled magma intrusion and released upward-migrated volatiles, in accordance with previous studies. Our inversion results bring new insights to the unrest. The intrusion, possibly in the form of a swarm of dykes or sills and dykes, arrived from depths below 8 km (b.s.l.) into the central part of the island just north of the Teide-Pico Viejo twin stratocones, causing a bulk density increase within this intruded volume related to a mass injection of roughly 2.3 · 10 11 kg. The intrusion propagated along the boundary between the basaltic core of the island, the Boca Tauce volcanic body and the more permeable and less compacted lower-density volcanic rocks (Figs. 7,8 and 9). The stalled intrusion released volatiles that disturbed the aquifer at the SW end of the caldera wall (injected mass of 2.9 · 10 10 kg). This unrest represents in our view a failed eruption, which is also in good agreement with the roughly 100 year repeatability interval of the recent historic eruptions. Our gravimetric model of the unrest indicates potential future magma ascent pathways and possible potential impact of the location of the assumed phonolytic magma chamber system.