New insights into La Pacana caldera inner structure based on a gravimetric study (central Andes, Chile)

La Pacana (central Andes, Northern Chile) is one of the largest resurgent calderas in the world, formed 4 Ma ago during an eruption with a VEI of 8.7. We undertake a gravimetric study to contribute new insights into the inner structure and evolution of this caldera. The La Pacana Bouguer residual anomaly is asymmetric and has an average amplitude of -12 to -14 mGal, which we interpret as being produced by the low-density intracaldera ignimbrite infill. A reinterpretation of the caldera stratigraphy plus the available geochronology suggests that the current shape of La Pacana was produced by the collapse of two nested calderas, roughly limited by the axis where the resurgent dome changes its orientation, with the oldest eruption in the southern part of La Pacana. The gravity data suggests that these southern and northern nested structures would have collapsed with downsag and trap-door geometries respectively, evidence of asymmetric subsidence. Intra caldera ignimbrite thicknesses were calculated with 2.5D forward models and show that the ignimbrite infill in its southern and northern parts reach ∼0.6-1.1 km and ∼2.5-3 km respectively. Using Gauss’s theorem, we calculate the volume of the intra caldera ignimbrite infill is ∼3,400-3,500 km3, in agreement with previous estimates and with models that show that the larger the caldera diameter, the larger the erupted volume.


Introduction
Resurgent calderas, also known as supervolcanoes, are the result of the largest volcanic eruptions on the geological record, which have been at least two orders of magnitude bigger than any eruption ever registered by humans (e.g., Newhall and Self, 1982;Mason et al., 2004;Miller and Wark, 2008). The large size of these eruptions (volumes larger than 450 km 3 and a Volcanic Explosivity Index (VEI) greater than 8) has attracted a growing global interest to understand the evolution and eruptive dynamics of the biggest calderas (Wilson, 2008), by studying their emitted products and the formation of their structure and geometry (Sparks et al., 2005;Acocella, 2007;Martí et al., 2008;Hardy, 2008). Due to the lack of observations of eruptions of this size, the caldera dynamics and the climatic impact of these supereruptions, including the effects over the human population, are poorly understood (e.g., Cole et al., 2005;Lipman, 1997;Self, 2006;Jones et al., 2007;Self and Blake, 2008). The largest ignimbrite f lareups in the world (e.g., Taupo volcanic zone, New Zealand; San Juan mountains, SW United States; Altiplano Puna Volcanic Complex APVC, southern central Andes of Bolivia, Chile and Argentina) are the places where the largest calderas in the world are located (Lipman, 2007;Rowland et al., 2010;De Silva, 1989a;De Silva et al., 2006;De Silva and Gosnold, 2007) and then constitute ideal study cases to investigate the dynamics of these supervolcanoes. The caldera inner structures are important in hazard assessment, as some of them show ground deformation likely related to magma intrusion (e.g., Wicks et al., 2006;Battaglia et al., 1999). Since caldera inner structures are normally covered by inf low ignimbrites, geophysical data are key elements to understand them and their dynamics.
In this work, we present a gravimetric study of the inner structure of La Pacana, the largest caldera in the Altiplano Puna Volcanic Complex (e.g., Gardeweg and Ramírez, 1987;De Silva, 1989a, b), which produced the fifth largest eruption ever registered in the geological record (e.g., Mason et al., 2004). Despite its size and importance, the existing geological studies (Gardeweg and Ramírez, 1987;Lindsay et al., 2001a, b) have several discrepancies regarding caldera limits, intra caldera ignimbrite thicknesses, collapse geometries and erupted volumes. We use gravimetry because it has been successfully used in caldera studies elsewhere to measure the subsurface caldera structure and deposit thickness -potentially resolving the discrepancies in previous interpretations (e.g., Rymer and Brown, 1986;Carle, 1988;Davy and Caldwell, 1998;Masturyono et al., 2001;Smith et al., 2006;DeNosaquo et al., 2009;Seebeck et al., 2009;Drenth et al., 2012). Based on both a gravity model and the reinterpretation of the ignimbrite stratigraphy, we suggest that the collapse geometry is significantly different than previously proposed, and that La Pacana caldera is made up of at least two nested calderas of different ages. We also provide estimates of the erupted volume and the collapse mechanism.

Geological Framework
The detailed geology of La Pacana caldera has been described elsewhere (Ramírez and Gardeweg, 1982;Gardeweg and Ramírez, 1985;Gardeweg and Ramírez, 1987;Lindsay, 1999;Lindsay et al., 2001a,b;Schmitt et al., 2002;Gardeweg and Lindsay, 2004). Here, we focus on the geological setting, the main structural elements including the resurgent dome and the caldera rim, and the caldera stratigraphy. La Pacana is the largest caldera of the Altiplano Puna Volcanic Complex of the central Andes (APVC), the largest Neogene ignimbrite f lare-up in the world, with an area of 30,000 km 2 (e.g., De Silva, 1989a, b;De Silva et al., 2006;De Silva and Gosnold, 2007). The APVC ignimbrite volumes are larger than 1,000 km 3 and their sources are several asymmetrical calderas that form huge depressions (e.g., De Silva and Gosnold, 2007), possibly with a tectonic control (Riller et al., 2001;. The APVC eruptive activity started 13 Ma ago (e.g., Schmitt et al., 2002) and occurred in short and intense eruptive pulses each ∼2-3 Ma (e.g., Salisbury et al., 2011). However, 2 Ma ago the volcanic activity decreased in magnitude and the eruptive style shifted to smaller and more effusive eruptions (e.g., De Silva and Gosnold, 2007;Salisbury et al., 2011). At a much larger, whole crustal scale, the APVC structure has been studied from a geophysical point of view by seismic tomography and receiver functions (e.g., Chmielowski et al., 1999;Zandt et al., 2003;Ward et al., 2014), electromagnetic methods (Schilling et al., 2006;Comeau et al., 2015), and gravity (e.g., Götze and Krause, 2002;Tassara et al., 2006;Riller et al., 2006;Del Potro et al., 2013). These studies show either the possible presence of magma body about 17 km below the surface (APMB, Altiplano Puna Magma Body) and an upper crust of felsic composition.
La Pacana caldera is located in the southern part of the APVC in the Antofagasta political region of northern Chile with its center at 67°25' W and 23°10' S ( Fig. 1 and 2) and is a 60 by 35 km asymmetrical caldera elongated in the NS direction. The caldera topographic margin is well exposed in its S part and has a maximum relief up to 1,000 m above the moat (Gardeweg and Ramírez, 1987). Gardeweg and Ramírez (1987) suggested that caldera subsidence was differential with the thickest intracaldera ignimbrite in the resurgent dome and decreasing towards both the north and south. However, the exact intracal-FIG. 1. Shaded SRTM topography of La Pacana caldera. The red and blue lines are the west and south caldera rims (Lindsay et al., 2001a). Open triangles show the caldera hanging wall, the orange line is Cordón La Pacana resurgent dome boundary, the red triangles are Miocene-Holocene volcanic centers (Ramírez and Gardeweg, 1982;Gardeweg and Ramírez, 1985;Lindsay et al., 2001a) and the black line is the Argentina-Bolivia-Chile border. Inset (colored from Lindsay et al., 2001a) shows the location of La Pacana caldera within the Altiplano Puna Volcanic Complex (APVC) in the central Andes.
dera ignimbrite thickness is unknown, as its base is not exposed. They proposed that the intra caldera ignimbrite maximum thickness is 0.9 km near the resurgent dome. Unlike the above mentioned study, the lack of a topographic margin in the N part of La Pacana led Lindsay et al. (2001a) to propose that the collapse mechanism was a trap-door geometry with a hinge region near the northern edge of the mapped caldera. These authors assumed a minimum thickness of intracaldera ignimbrite of 2 km for La Pacana as a whole, which would be thicker in the southern part of the caldera, an infill discrepancy of more than  (Lindsay et al., 2001a;Gardeweg and Lindsay, 2004). The rhyolitic centers within the southern part of the caldera moat (red bodies south of Cordón La Pacana) have been modified from the original map. Dashed blue lines and text labeled as NW, NE, SE and SW are the trace of the modeled cross sections in figures 4 and 5, and the dashed red line shows the inferred limit between the southern and northern caldera collapses (see discussion section).
For the most part, La Pacana overlies a sequence of late Miocene pre caldera ignimbrites and volcanic rocks, and locally, upper Paleozoic volcanic and lower Paleozoic sedimentary basement (Ramírez and Gardeweg, 1982;Gardeweg and Ramírez, 1985;Gardeweg and Ramírez, 1987;Bahlburg et al., 1988;Breitkreuz and Zeil, 1994;Lindsay et al., 2001a) (Fig. 2). Due to the complexity of the caldera and the lack of dissection, the La Pacana collapse has been proposed to be produced by the eruptions of either Atana (4.11-4 Ma) (Gardeweg and Ramírez, 1987), Atana and Toconao (4.65 Ma) (Lindsay et al., 2001a,b;Schmitt et al., 2002), Pujsa (5.7 Ma) and Atana (Pavez et al., 2008) ignimbrites, or all of them (De Silva and Gosnold, 2007) (Fig. 2). Regardless of the erupted ignimbrites, the eruptions resulted in ponding of these units within the collapse area into a thick sequence of intracaldera ignimbrite and a smaller outflow volume. Some of the ignimbrites inside the caldera might have been erupted from other calderas in the area (Lindsay et al., 2001a). Post-eruption resurgence of La Pacana led to the formation of an uplifted asymmetric elongated dome of intracaldera tuff known as Cordón La Pacana. This dome occupies 24% of the area within the topographic depression with long and short axis of maximum lengths of 48.5 and 12 km trending in the NW-SE and NE-SW directions respectively ( Fig. 1), with its highest part reaching ∼5,200 m above sea level, almost 1 km above the caldera moat lying at ∼4,200-4,500 m. The dome is broken by several NW-SE and SW-NE striking normal faults produced by its resurgence (Gardeweg and Ramírez, 1987). This asymmetrical geometry and a reinterpretation of the domes that intrude the southern part of the resurgent dome (Arenoso, Chivato Muerto, and Chamaca domes, red bodies in the caldera moat in figure 2, syn-caldera domes hereafter) and their available K-Ar ages (Gardeweg and Ramírez, 1987;Lindsay et al., 2001a) led Pavez et al., (2008) to propose that La Pacana could be at least two nested calderas (the inferred limit is the dashed red line in figure 2) produced by two major eruptive events (see discussion for further details). Finally, post caldera volcanism lasted from 4 up to 2 Ma (Gardeweg and Ramírez, 1987;Lindsay et al., 2001a) and occurred as dacitic dome and stratovolcanoes extrusions along the western, southern, and eastern edges of the resurgent dome (Lindsay, 1999).

Data Acquisition and Processing
71 gravity stations were acquired in a survey in January-February 2009 with a LaCoste & Romberg G411 gravity meter and a pair of Leica SR20 single frequency geodetic receivers with Leica AT501 antennas. Due to the large distances in the survey area, the gravity base station was measured at the beginning and the end of the daily surveys and several cross-loops were measured using relative stations tied to the reference benchmark. Additional observations were provided by 20 stations acquired by Pavez (2005) and 647 gravity stations acquired between 1,982 and 1,994 from the central Andes gravity data base . Due to its pre GPS date, the elevation benchmark of each station of the latter database had to be validated before it was merged with the newest data.
The dataset was reduced using standard procedures (Telford et al., 1990;Blakely, 1996). The earth tide was corrected with Longman's formula (1959), while the drift correction was calculated with the CGxTOOL software (Gabalda et al., 2003). Resulting drifts were lower than 0.3 mGal/day, in agreement with literature (Blakely, 1996). The 2009 data set was tied to the central Andes database, which is part of the IGSN71 gravity net (Morelli et al., 1974). On the other hand, GPS benchmarks were tied to the Caltech Chajnantor (CJNT) continuous geodetic station (Simons et al., 2010). The GPS bases for each day of the survey were then processed in the Spectra Precision Survey Office software in static mode tied to CJNT, decimating the data to the common sampling rate of 0.033 Hz. Then, the GPS rover of each gravity station was processed in kinematic or fast static mode tied to its daily GPS base. Baselines were shorter than 25 km and measurement times of each GPS station were 7 minutes on average, yielding a vertical accuracy of less than 1 m for the 2009 data. As an example, for 3 and 25 km GPS base lines, horizontal accuracy was less than 0.3 and 0.8 m respectively. Due to the lack of dual frequency GPS positioning of the central Andes gravity data base, the stations where repositioned with pre existent kinematic GPS tracks and their elevations were recalculated using the 3 arcsec Shuttle Radar Topographic Mission (SRTM) digital elevation model (DEM) (Farr et al., 2007). This DEM has a spatial resolution of 90 m/pixel and was leveled by adding the difference between this data set and the elevation of all the geodetic GPS stations and tracks (35 m), and then used to extract the elevation benchmarks of the latter gravity database.
The accuracy of the different gravity data sets is variable. The Lascar volcano gravity data (Pavez, 2005) accuracy is better than 0.05 mGal, while the accuracy of the 2009 survey varies between 0.05 and 0.1 mGal. The accuracy of the elevation benchmarks of the gravity stations can be compared converting the former figures to the free air gradient of 0.3086 mGal/m. This means that the elevation benchmarks of the 2009 data are better than 0.25 mGal, which are the expected uncertainties in a regional survey like this. Hence the combined 2002-2009 data accuracy is better than 0.22 mGal. Schmidt and Götze (2006) data are not as accurate as the former because the difference between the GPS and SRTM elevation benchmarks has a standard deviation of 13 m, therefore the gravity related error is 4 mGal of free air correction, still an acceptable figure for this regional work.
The Bouguer residual anomaly separates anomalies produced by bodies not accounted for by a simple density reference model of the crust and was calculated with standard expressions (Moritz, 1984;LaFehr, 1991a;Blakely, 1996), implemented in the Geosoft Oasis montaj software. The topographic correction was calculated with the SRTM DEM using standard algorithms (Nagy, 1966;Hammer, 1939;Kane, 1962) and the Bouguer correction included the Bullard B correction (LaFehr, 1991b) to account for the Earth's curvature. The terrain density for the Bouguer correction is one of the most critical parameters in gravity reduction because the resulting anomalies can be over or underrepresented depending on the density contrast, especially when the lithological variations are large such as the study areas with multiple ignimbrites and facies. On the other hand, the nonuniqueness separation of the regional/residual field associated to long-wavelength-deep-seated-bodies/ shorth-wavelength-shallow-geological-features, can also lead to over or underrepresentation of geological anomalies (e.g., Gupta and Ramani, 1980;Simpson et al., 1986;Blakely, 1996). A simple form to calculate a reasonable density and isostatic model is to minimize the Bouguer residual anomaly by linear least squares, to reduce its correlation with the topography (Nettleton, 1939;Parasnis, 1962;Torge, 1989) a method used in volcanic studies elsewhere (e.g., Tiede et al., 2005;Pavez, 2005). The inverted density is 2.2 g/cm 3 , which is a reasonable value for a lowdensity ignimbrite covered area. The resulting Bouguer residual anomaly was finally gridded with the minimum curvature algorithm (Briggs, 1974) with a cell size of 3 km.

Bouguer residual anomaly
The Bouguer residual anomaly shows a strong and asymmetrical NS elongated negative anomaly (shades of green and blue in figure 3) which is broader than La Pacana caldera boundaries, with an average amplitude of -12 to -14 mGal reaching -21 and -40 mGal in the southern and northern parts of the caldera respectively, the global minimum located in the northern part of the resurgent dome (Fig.  3). This signal is surrounded by strong positive anomalies with maximum amplitudes of 30 and 45 mGal west and east of La Pacana. All of these signals have amplitudes at least 3 to 9 times above the error level. Different interpretations on the origin of the low density-negative gravity anomaly have been proposed including the APMB or a batholith (De Silva et al., 2006), brittle felsic upper crustal rocks  or low-density ascending diapiric bodies rooted in the APMB in the northern part of the APVC (del Potro et al., 2013). However, as the observed negative anomalies are spatially related with La Pacana caldera, we consider that the most likely source of this signal is the low-density intracaldera ignimbrite facies, in agreement with studies in caldera environments elsewhere (e.g., Rymer and Brown, 1986;Carle, 1988;Masturyono et al., 2001). On the other hand, relative positive anomalies such as the signal in the middle of the caldera and spatially related with volcanic edifices (Fig. 3) are interpreted as being produced by dense intrusive bodies (e.g., Rymer and Brown, 1986;Malengreau et al, 1999;Gudmundsson and Högnadóttir, 2007), due to the crystallization of the densest mineralogical phases during magma ascent and emplacement. Finally, positive anomalies outside of the caldera are interpreted as the response to the uplift of dense Paleozoic basement.

2.5D Forward Modeling
In order to quantify the Bouguer residual anomaly interpretation, we calculated 2.5D forward models along two transects. These sections allow an estimate of the intracaldera infill thickness, margin geometry, and geometry of the underlying bodies. These profiles (black lines in figure 3) were chosen to traverse gravity anomalies of geological interest in places where the contribution to the gravity field is mostly associated with the newer gravity data, ensuring that uncertainties are below 1 mGal. The forward   (Lindsay et al., 2001a). Open triangles show the caldera hanging wall. The dashed brown line is the inferred border between the N and S parts of La Pacana. Black, white and grey circles are gravity stations acquired in the 2009 summer, by Pavez (2005) and Schmidt and Götze (2006). Red triangles are Miocene-Holocene volcanic centers (Ramírez and Gardeweg, 1982;Gardeweg and Ramírez, 1985;Lindsay et al., 2001a). White dashed lines delimit positive anomalies interpreted as being produced by Paleozoic basement. White box is an example of a relative positive anomaly interpreted as being produced by intrusive bodies beneath the volcanic edifices. Black lines indicate the selected profiles where gravity forward modeling was carried out. White letters indicate the directions of the forward models in figures 4 and 5. The grey box is the extent of the grid used in the calculation of the intra caldera mass deficit.
modeling was carried out in the software ENCOM ModelVision Pro, based on the Talwani's algorithm implemented by Cady (1980) for the case of 2.5D bodies of limited extension and symmetric in the third coordinate and with different densities. The densities of the modeled bodies were measured from field samples, and in the case of intrusive bodies, extrapolated from dome densities ( Table 1). The polygon shapes were constrained with geological information when available, such as mapped caldera structures or borehole data. In order to reduce the strong non-uniqueness of the forward models (e.g., Roy, 1962), intrusive bodies shapes are kept as simple as possible, thus simple shapes such as ellipsoids approximated some of them. Therefore, we consider the modeled geometries a simplified model of the real caldera structure. Based on the arguments of the previous section, we make the assumption that the negative gravity signal is produced by low-density intracaldera infill and that the underlying basement has no density contrast with respect to the density reference model. Finally, we state that there are minor differences of 100-200 m in the thickness of the modeled bodies at the intersection points of the profiles, which are produced by considering the 3D caldera structure as 2.5D, which do not change the overall interpretations of this work.

SW-NE Cross Section
This section (Fig. 4) traverses the southern part of La Pacana caldera, where the gravity signal reaches up to -21 mGal, thus the deepest area of this section of La Pacana. The gravity anomalies are always negative, but reach relative maximums within the caldera moat (-2 mGal) and in the post caldera volcanoes (-7 mGal), while the minimum values are below the resurgent dome. The profile is modeled with five bodies, two with positive densities interpreted as shallow intrusive bodies feeding domes and stratovolcanoes, the third with a negative density contrast interpreted as La Pacana intracaldera ignimbrite infill, and two with negative density contrast required to fit the gravity signal outside of the caldera and interpreted as undifferentiated low density ignimbrites. The intracaldera thicknesses are variable, with a minimum of 0,6 km below the caldera moat and a maximum of 1,1 km in its deepest part below the resurgent dome with some f luctuations in between (Fig. 4). The possible collapse geometries (Lipman, 1997;Cole et al., 2005) arising from the modeled ignimbrite infill are described in the discussion section. We speculate that the intracaldera ignimbrite infill in the southern part of the caldera is made up of several ignimbrites besides Atana (see discussion section), but due to the lack of density contrast we can not determine the boundary between the different volcanic deposits in the gravity model. The ignimbrite infill is modeled with respect to the described topographic margin, limited by normal faults according to analog models (Acocella, 2007;Hardy, 2008;Martí et al., 2008) and with their sense of slip inferred from these models (Fig. 4), therefore more complex structures such as the caldera collapse collar (Lipman, 1997) are not taken into account Although andersonian normal faults are expected to dip ∼60º, the modeled dips are shallower and steeper along the interpreted structure. Finally, the forward model shows that the roof of shallow intrusive bodies is located at ∼1.5-2 km below the surface and the ignimbrite sheets. These intrusions are broader in wavelength than the fed syn-caldera domes and post caldera stratovolcanoes, and given their densities they are most likely of dacitic composition. Due to the strong non-uniqueness in gravity forward models and the lack of further constraining data, we have assumed that the material between the domes and the intracaldera infill does not have a density contrast with respect to the density used in the Bouguer residual anomaly calculation.

NW-SE Cross Section
The second cross section (Fig. 5) was chosen to traverse the northern part of La Pacana through the caldera moat and the resurgent dome. The gravity anomaly decreases from 18 mGal in the SE part, to -30 mGal in the NW, with some fluctuations in between. Three bodies were considered to fit the gravity signal, which account for the intracaldera infill, an intrusive beneath the resurgent dome and the pre-volcanic basement. Unlike the previous section, we model part of the pre-volcanic basement as it outcrops east of La Pacana (Gardeweg and Ramírez, 1985) and the gravity signal of the uplifted basement is required to fit the data. The modeled ignimbrite thickness increases towards the NW reaching a maximum of ∼2.5-3 km, but the NW limit of the ignimbrite infill cannot be properly determined with the limited data set. The intrusive body below the resurgent dome produces the local maximum in the gravity signal as a result of being located close to the highest part of the resurgent dome; hence its

322
New iNsights iNto La PacaNa caLdera iNNer structure based oN a gravimetric study (ceNtraL aNdes, chiLe) emplacement could have been responsible for the dome uplift as shown by overpressure analogue experiments (Acocella, 2007).

La Pacana caldera depth infill estimates and limits: evidence for the presence of nested calderas
The modeled intracaldera infill thickness is validated with borehole data in the southern part of the caldera moat (M. Gardeweg, personal communication, 2009) on which the pre volcanic basement lies 0.6 km below the surface. Thus, the ignimbrite infill in this area has a maximum depth of 1.1 km in the resurgent dome and decreases towards the southern part of La Pacana (Fig. 4), while the ignimbrite infill reaches 2.5-3 km in its NW part (Fig. 5). The former ignimbrite model is more akin to the 1 km of intracaldera infill that Gardeweg and Ramírez (1987) proposed for the central part of the caldera than the 2 km assumed by Lindsay et al. (2001a). The caldera northern boundary cannot be properly delimited with the available data set as the large scale negative gravity signal (Fig. 3) does not show clear boundaries and is not well constrained by the newer gravity data set.
A reinterpretation of the La Pacana available stratigraphy and geochronology lead us to consider a different origin for the modeled intracaldera ignimbrites than the inferred interpretations arising from Gardeweg and Ramírez (1987) and Lindsay et al. (2001a). The southern part of Cordón La Pacana dome is intruded by three syn-caldera domes (red units within the caldera moat in figure 2), for which only a single age K-Ar age of 4.8+0.2 Ma is available Ramírez, 1985, 1987) that lead these authors to suggest that these units are pre-caldera collapse as they are older than the Atana ignimbrite. However, Lindsay et al. (2001a) pointed out that as they lie over the southern part of the resurgent dome, they must be younger than the latter unit. The age of these syn-caldera domes overlaps with the only Atana ignimbrite K-Ar age for the southern part of the resurgent dome of 4.5+0.4 Ma (Gardeweg and Ramírez, 1987). But most of the available Atana K-Ar ages are younger and span ranges between 4.2 and 3.7 Ma with an average of 4.0 Ma (Gardeweg and Ramírez, 1987;De Silva, 1989b;Lindsay et al., 2001a), while its U-Pb ages are 4.11+0.2 Ma (Schmitt et al., 2002) therefore we consider the 3.7-4.2 Ma time span as the most representative Atana age. This implies either the syn-caldera domes were intruded into the resurgent dome quickly after the emplacement of the latter or as they lie within the caldera moat, they are older than the resurgent dome in the southern part of La Pacana caldera. As the accepted ages for the eruption of the Atana ignimbrite are younger than the syn-caldera domes, we suggest that the southern part of La Pacana was produced by an older eruption and collapse than the Atana eruption that generated its northern part as older volcanic edifices do not intrude the resurgent dome in this area (Ramírez and Gardeweg, 1982;Gardeweg and Ramírez, 1985), Hence, the observed caldera geometry would be the result of a nested collapse, whose limit is inferred to be roughly the zone where the resurgent dome changes its orientation from NW-SE to NE-SW (brown dashed line in figures 2 and 3). Therefore most of the intracaldera infill modeled in figure 4 is made up of Atana and older ignimbrites, but we cannot assess with the available data whether any of the pre-Atana ignimbrites like Toconao (Fig. 2) were erupted from this area. We acknowledge that this interpretation is highly speculative as it is based in single K-Ar ages for each of the geological units and the low spatial resolution of the gravity data set. The lack of density contrast between the ignimbrites hampers the ability to calculate gravity gradients that could be used to determine density boundaries between the two nested calderas. In the absence of more observations and data, we have not attempted to further redefine the caldera limits besides the first order observation that the boundary may be where the resurgent dome changes strike. On the other hand, ellipitical calderas have been proposed to be elongated parallel to the minimum horizontal stress (Holohan et al., 2005), in agreement with regional tectonics that show an horizontal σ 1 of west-east direction during the Pliocene (González et al., 2009). However, La Pacana elliptical shape could also be produced by the nested collapse of the smaller calderas, which hampers the assessment of whether the caldera set final shape is tectonically controlled or not.

Collapse mechanisms and geometry
The different gravity signals inside La Pacana caldera provide new insights about the caldera collapse mechanism and geometry. The superposition of the intrusive gravity signals inside La Pacana complicates the interpretation of the collapse geometry as the Bouguer residual anomaly is produced by multiply bodies; hence our data can not be used to support complex collapse geometries such as a piece-meal model in which several blocks fall asymmetrically during the caldera collapse. Gardeweg and Ramírez (1987) suggested a central ring fault around the resurgent dome, while Lindsay et al. (2001a) noted the lack of a northern topographic border and proposed a trap-door geometry with a northern hinge zone. Based on the previous section, we assume that the caldera is made up of two different collapses, hereafter PS (Pacana south) and PN (Pacana north). The presence of a negative anomaly close to PS central part (Fig. 3) suggests a downsag collapse geometry (Cole et al., 2005), with maximum subsidence near the caldera center. Compared with the earlier Gardeweg and Ramírez (1987) and Lindsay et al. (2001a) proposals, our model for PS caldera is more consistent with the former. On the other hand, as the forward model (Fig. 4) shows that the gravity signal is deeper towards PN northwest zone, we suggest a trap-door mechanism for this caldera with a hinge zone in the area where the resurgent dome changes its orientation (brown line in figure 3).

Mass deficit and erupted volume (VEI)
One of the key parameters in all caldera studies is the volume of the ignimbrite infill and outf low, which allows to estimate the VEI of the associated eruptions (e.g., Miller and Wark, 2008). The volume estimation of old eruptions as in the APVC is a nontrivial task due to the lack of dissection that hampers estimation of the ignimbrites' thicknesses and the petrological similarity of ignimbrites erupted from different calderas. This calculation is based on simple assumptions such that the intracaldera and outf low volume deposits are approximately equal; therefore the volumes are only order of magnitude estimates (e.g., Mason et al., 2004). These assumptions do not hold for the APVC calderas (Salisbury et al., 2011). For La Pacana caldera, the previous work has presented different erupted volumes, based on the assumption of different geometries for the caldera inner structure (Lipman, 1997). Gardeweg and Ramírez (1987) estimated the erupted volume to be ∼900 km 3 , while Lindsay et al. (2001a) estimated ∼2,500 km 3 . Unlike the previous studies, we propose to use Gauss's theorem to calculate the mass deficit produced by the caldera eruptions and the associated volume, which has been used in other calderas (Campos-Enríquez et al., 2005). According to Gauss's theorem (e.g., Hammer, 1945;LaFehr, 1965), the mass M related to a gravity anomaly Δg measured over the Earth's surface (S) is with G the gravitational constant (6.67*10 -11 Nm 2 /kg). This formula is discretized for cells with constant area S as where g i is the gravity anomaly of the i grid cell. The main advantage of Gauss's theorem is that the mass estimate is independent of the anomalous body geometry. Therefore, any VEI calculated with the combined mass of the coeruptive ignimbrites must be regarded as a total value that cannot discriminate different eruptions. However, it posses some disadvantages, such as a proper regional, residual gravity fields separation, which could over or underrepresent the mass related to the gravity anomalies. Another disadvantage is that the exact surface that contains La Pacana caldera is not regular, so if the calculation is carried out with a square grid, the resulting mass deficit figure has to be reduced to take into account an estimate of the percentage of cells that do not contain the caldera gravity signal. The only way to improve this procedure (and beyond the scope of this paper) is to calculate a 3D model of the caldera structure and estimate the mass directly from the volume of the modeled polyhedrons.
The gravity grid used to calculate the mass deficit is shown within the grey box in figure 3. As part of the gravity signal of the ignimbrites is outside the caldera, the box extends some distance from the caldera topographic margin. The calculated mass deficit of La Pacana is -1.08*10 15 kg, but as the grid is bigger than the caldera extent, the previous figure has to be reduced 20%, with a resulting mass deficit of -8.6*10 14 kg. We calculate the volume with the same density contrast of -300 kg/m 3 that was used in the caldera intracaldera forward models, with a resulting ignimbrite infill of 2,868 km 3 . The outflow volume has been estimated between 240 (Gardeweg and Ramírez, 1987) and 270 km 3 (Lindsay et al., 2001a), so we assume an average of 255 km 3 . An additional outf low is recorded below the Salar de Atacama basin (Jordan et al., 2002), with average thicknesses of 0.2 km over an area of 3,000 km 2 (Jordan et al., 2002) that yields a volume of 600 km 3 . Nonetheless, a conservative estimate is that the ignimbrites are present below 1/2 and 2/3 of the Salar de Atacama basin, thus reducing the volume to ∼300-400 km 3 . Therefore, the caldera outflow is ∼555-655 km 3 . The final volume of the caldera infill and outflow resulting from these estimates is ∼3,400-3,500 km 3 , larger than by Lindsay et al. (2001a) and de Silva and Gosnold (2007) for the combined Pujsa, Toconao and Atana ignimbrites (∼2,800 km 3 ). Thus, the resulting volume of the PC and PS calderas makes them some of the largest in the world, comparable to Yellowstone, Toba and La Garita (Mason et al., 2004). Assuming a vesicularity of 35% (Lindsay et al., 2001a), the ignimbrites DRE (dense rock equivalent) is ∼2,200-2,300 km 3 .
The associated VEI is calculated with a modification of the the original scale by Newhall and Self (1982), VEI = log 10 (m)-7 (Pyle, 1995), with m the mass in kg of the caldera ignimbrites. The mass of the outflow and the caldera infill is ∼1.02-1.05*10 15 , thus yielding an average VEI for both ignimbrites of ∼8, a figure slightly smaller than the Mason et al. (2004) calculation for the Atana eruption, but anyway one of the largest in the world. We note that Mason et al. (2004) interpreted their VEI as being produced by the single eruption of Atana ignimbrite, while our calculation considers all the erupted products regardless of the eruption. This discrepancy can also be explained by the smaller density values used in this work (1,900 kg/m 3 ) while Mason et al. (2004) used larger values (2,200 -2,400 kg/m 3 ), thus yielding different estimates. Regardless of the erupted mass and VEI the caldera set considered as a single body (as we cannot discriminate the volumes of La Pacana southern and northern parts) follows the linear logarithmic trend that the larger the caldera radius, the larger the mass deficiency (mass of the erupted volume) (Fig. 6). FIG. 6. Comparison of caldera diameter and mass deficiency. Blue circles are data from Campos-Enríquez et al. (2005) and Yokoyama (1987). Grey circles are caldera erupted volumes converted to mass def iciency with a density contrast of 300 kg/m 3 (Acocella, 2007). Red and yellow circles are La Pacana caldera with the mass deficiency calculated by means of Gauss's theorem and from Mason et al. (2004). La Pacana radius was calculated as the mean of the north-south and east-west diameters (Gardeweg and Ramírez, 1987).

Conclusions
Despite the low spatial resolution of the gravity data set, first order inferences can be made about the inner structure of the La Pacana caldera that can be compared to different geological interpretations based on previous work. Our reinterpretation of the caldera stratigraphy led us to speculate that its southern part was produced by a different eruption and collapse than its northern zone, therefore La Pacana might be made up of two nested smaller calderas. The modeled maximum intracaldera ignimbrite thicknesses are 0.6-1.1 km and 2.5-3 km for the northern and southern parts that collapsed with downsag and trap-door geometries, different than previous models but more consistent with Gardeweg and Ramírez (1987) than Lindsay et al. (2001a) proposals. The average erupted volume of La Pacana is ∼3,400-3,500 km 3 (VEI=8) for the combined eruptions, in agreement with models that show that the larger the caldera area, the larger the erupted volume. Further work that will help to better understand the caldera dynamics and clarify its evolution include detailed geochronological work, and denser geophysical surveys that incorporate gravity, aeromagnetic and magnetotelluric data.