The PCCVC (40.5°S, 72.2°W) is a cluster of Pleistocene-to-Holocene volcanic vents aligned in NW–SE trend oblique to the main volcanic front of the Chilean Andes Southern Volcanic Zone (SVZ; Fig. 1a) (Lara et al. 2006a). Uniquely, the PCCVC erupts high-MgO basalt together with rhyodacite and rhyolite products, unlike those observed in most basaltic and basaltic andesite SVZ volcanoes (López-Escobar et al. 1981; Hickey et al. 1986; Gerlach 1988; Tormey et al. 1991). The products consist of lava flows, lava domes, pumice fall deposits and lahars sourced from at least 27 vents from the Late Pleistocene (16.5 ka) to historical times (Katsui and Katz 1967; Lara et al. 2004, 2006a; Lara and Moreno 2006; Singer et al. 2008).

At least six Holocene tephra fall deposits have been identified from the PCCVC, with bulk volumes ranging from 0.1 to 4.32 km³ and a magnitude of about 5-6 (VEI 5; Naranjo et al. 2017; Toloza et al. 2020). A series of 0.6-0.7 ka BP sub-Plinian phreatomagmatic eruptions destroyed the summit of Puyehue volcano forming the Mil Hojas sequence, totaling 4.3 km³ bulk volume (Naranjo et al. 2017; Alloway et al. 2021). In 1921-1922 a sub-Plinian eruption was followed by an effusive stage and erupted ca. 0.3 km³ of DRE volume (Steffen, 1922; Krumm 1923; Katsui and Katz, 1967; Singer et al., 2008). Another sub-Plinian eruption was triggered 38 h after the May 22, 1960, M_w 9.5 earthquake (Veyl, 1960; Katsui and Katz, 1967; Lara et al. 2004, 2006b; Lara and Moreno, 2006; Singer et al. 2008), erupting ca. 0.57 km³ tephra (DRE ca. 0.25 km³; Chapron et al., 2006; Lara et al. 2006a, 2006b). Both eruptions are characterized by rhyolitic to rhyodacitic bulk compositions (68.90-70.01 wt. % SiO₂). The tephra deposits from the two sub-Plinian eruptions are clearly recognizable in outcrops along the 215 international highway as two discrete, decimeter-scale layers containing partially weathered pumice lapilli. Petrological evidence suggests that the historical eruptions have all been fed by the same long-lived mush system (Seropian et al. 2021).

Seismic swarms from a volcano-tectonic source occurred in May 1994, March 1999, and May 2007 (Lara et al. 2006a) and two uplift events occurred within the PCCVC between 2008 and 2010 (Jay et al. 2014). The 2011-2012 explosive eruption began at a site known as “Las Sopas” that formerly produced sulfur emissions, hot springs, steaming-ground, and acid-sulfate alterations (Sepúlveda et al. 2005).

The June 4, 2011, PCCVC eruption was preceded by three ground uplift periods between 2003 and 2010 (Delgado, 2021) and a seismic unrest that began in December 2010 (Basualto et al. 2023). It consisted of an increase in the number and magnitude (M > 3) of volcano-tectonic, long-period, and hybrid seismic events that became shallower (1-4 km depth) until April (Bertin et al. 2015). A source of deformation was observed between March and May 2011, covering 40 km² and accounting for a 6 cm uplift (Jay, 2014). In May 2011, the magnitude of earthquakes reached M~4 and one day prior to the eruption their frequency was 25 events per hour (Lara et al. 2012).

The eruption started around 18:45 UTC on June 4, 2011, and initially produced a sustained sub-Plinian eruption column (~14 km above sea level, asl). The initial phase lasted 27 hours and the tephra fall caused severe disruption over the Argentinean Patagonia (Collini et al. 2013; Schipper et al. 2013; Jay et al. 2014; Murriello and García, 2023) and around PCCV in Chile (e.g., Marin et al., 2020; Vergara-Pinto and Marin, 2023; Vergara-Pinto and Romero, 2023). In the following days, the eruption column oscillated between heights of 3–10 km asl fed by the waning explosive phase (Castro et al. 2013). Lava effusion started on June 15 (Bertin et al. 2015), which was coincident with explosive activity characterized by mixed gas and ash jetting punctuated by Vulcanian blasts with eruption columns up to 6 km high (Schipper et al. 2013). During the eruption, the explosivity correlated directly with the number of earthquakes per hour. Most of the tephra dispersal occurred towards the NNE and ESE as a result of the prevailing winds. The effusive phase had discharge rates up to 70 m³s^-1 and finished by mid-2012 (Bertin et al. 2015). The rhyodacitic lava flow covered ~7 km² with the extent limited by a graben and the 1960 lava flow to the south. While the 1960 eruption formed a 4 km-length fissure system, the 2011-2012 eruptions were centralized, forming a 0.5 x 0.9 km-diameter pyroclastic cone. Northeast of this cone, proximal tephra fall sequences, and composite PDC-tephra fall-ballistic deposits reached ~3 m and > 15 m-thick, respectively. Since the eruption, at least three ground uplift events have occurred between 2012 and 2019 at rates from 6 to 45 cm/yr, resulting in >1 m total uplift (Euillades et al., 2017; Delgado, 2021).

The 2011- 2012 tephra fall deposit was characterized in detail by Pistolesi et al. (2015). The pyroclastic units from the main sub-Plinian phase and the activity up to June 11 correspond to Unit I (layers A to F, deposited by 4-5 June), Unit II (layers G to H, deposited by 5- 6 June) and Unit III (layers K1 to K5, deposited by 7-11 June). Within 25 km from the vent, the sequence is dominated by multiple lapilli-bearing fall layers (Units I and II) overlain by cm-thick, light gray, ash-rich fall layers (Unit III). In the most proximal outcrops, along the slopes of the PCCVC, a fourth and younger unit can also be recognized and is represented by very fine-grained, thin (few millimeters to 1-2 cm) white ash related to the eruptive activity continuing to the end of the eruption in 2012.

As shown in Table 1, several authors have estimated the bulk and DRE volume of the tephra fall deposit.

The most proximal locations observed E and NE from the vent (~1.2 km) consist of a coarse grained, lapilli-to-bomb sized, pumice tephra fall deposit of ~3 to 3.2 m in overall thickness. These deposits are well bedded and have alternate finer and coarser pumice fall and ash layers. Complex sequences lying NE from the vent reflect the interplay between the coincident tephra fall and pyroclastic density currents during the sustained phase of the eruption.

At medial downwind locations, the tephra fall deposit consists of 3 to 9 pumice lapilli-to-ash layers (Fig. 1). From base to top, the thickest described section (Fig. 1a; 215 highway, 57 cm-thick) starts with layer A consisting of ~7 cm-thick ungraded fine white sub-rounded lapilli pumice; largest-sized pumice fragments are highly vesicular and up to 0.7 cm in diameter. Above, layer B is a ~14 cm-thick non-graded fine-to-medium white pinkish lapilli pumice layer with largest clasts up to 2.5 cm in diameter. Layers A and B contain 5-15 % angular black cognate lithics (0.1-0.3 cm diameter). Layers C and D (8 and 5 cm-thick, respectively) are inversely graded and contain fine-to-medium clear-brown to yellowish blocky lapilli pumice. To the top, the sequence has two layers, named E and F, with thicknesses of 4 and 3 cm, respectively; these layers are composed of ungraded, medium yellowish subrounded or blocky lapilli with the largest particles reaching 1.5 and 2.0 cm in diameter. Increasing concentrations of angular black lithics are observed towards the top, some up to 1 cm in diameter. While layer H consists of ~2 cm-thick coarse ash and fine lapilli with a white-grayish color, the package of K layers is ~4 cm-thick and contains fine yellowish lapilli pumice with overall reverse grading. The uppermost layer (L) consists of ~2 cm of finely laminated gray ash. According to the interpretation of Pistolesi et al (2015), layers A-F comprise Unit I, while Unit II are layers G and H. All the layers in K form Unit III, whereas Unit IV corresponds to layer L. The integrity of the so-called Unit I is generally well-preserved in the farthest exposures of the medial sites (Fig. 1b, c), however, layers A and D-F thin out rapidly and the contact between B and C becomes diffuse.

The medial crosswind deposits, north of Puyehue-Cordón Caulle are different from the downwind exposures: the basal layer is a massive coarse brown pumiceous ash (A-F; ~1.5 cm-thick). On top, layer G (~3.0 cm-thick) is composed of reverse-graded fine-to-medium grained white lapilli pumiceous fragments (~ 50%), black volcanic glass (~ 40 %), and dark angular lithic fragments (~10 %) from 0.4 to 1.0 cm diameter. However, this layer has only been found up to 6.5 km N from the vent (Fig. 1d). Above this layer, the Unit III sequence consists of at least five layers (totaling >5 cm), characterized by alternating white fine ash and dark-black glassy angular fragments. To the lower part of the sequence, the contacts are sharp, while towards the top they are mainly gradational.

At distal locations, the downwind tephra fall consists of a single homogeneous white or grayish ash fall layer with some color grading in the middle of the vertical exposure (Fig. 1e, f) with the thickness varying from a few mm to a few cm.

The integrity of the proximal-to-medial sections agrees with the stratigraphy described by Pistolesi et al (2015). Moreover, we did not notice evidence of remobilization or deposit transformation at selected sites for this study. Thus, we suggest all our observations well represent the primary tephra fall deposit and are suitable for volume estimates.

The full-dataset isopach map (Fig. 3) shows E-SE distribution of the tephra deposit, mainly over Patagonia Argentina. During the eruption, most of the Argentinean Pampa and some populated cities such as Neuquén were affected by the distal deposition of ashfall (Fig. 3a). The thinnest measured isopach was 1 mm thick and encompasses nearly 70,000 km². Ash traces, fine enough to be transported over longer distances may have fallen in the Atlantic Ocean, thus accounting for a volume difficult to constrain. At intermediate locations (up to ~230 km downwind), the tephra deposit reaches between 2 and 6 cm-thick to the E-SE. That portion of the deposit occupies an area of ~20,000 km² and its shape is highly irregular; it develops two main dispersal axes to the NE and E-SE, because of complex high-altitude and surface wind regimes and/or the local influence of topographic features. The thickest isopachs (> 6 cm) display regular, ellipsoidal shapes and are distributed obliquely to the 215 and 231 international highways, over the town of Villa La Angostura (Fig. 3b), up to 40 km in the downwind direction. We adopt the nomenclature of Pistolesi et al. (2015) to describe the internal deposit architecture.

In all the isopach map iterations, increasing the FCP number in the different datasets yields an increasing number and complexity of isopachs. Isopach maps constructed from only 10 FCPs produces 3 or 4 rough isopachs, in contrast to ≥10 isopachs being drawn when the datasets include ≥110 FCPs (Supplementary Table S4). Increasing the number of FCP points used corresponds to isopachs changing from elliptical to irregularly shaped contour lines (Fig. 4a-e) which are considered more realistic presentations of the natural deposition of the tephra. For example, the isopach map based on the 150 FCP dataset is fairly like that shown in Fig. 2a, which contains the full FCP dataset. On the other hand, the square root area v/s thickness regression illustrates higher variability for smaller FCP datasets, while larger FCP number datasets generally display a more conservative distribution (Fig. 4f). Considering the full dataset as the most realistic (or exact isopach areas), the absolute percent error of the isopach square root areas obtained with smaller datasets range between 1 and 63%, with average errors in the range 20-33 % and mostly concentrated in the distal isopachs (<0.01 m-thick; Fig. 4).

The bulk volume estimates of the full dataset isopach maps are 1.1 and 1.2 km³ using the Weibull and exponential methods, respectively (Supplementary Table S3). Therefore, our best volume estimate using both methods is 1.15 ± 0.05 km³. However, important volume differences are found using the variable FCP number dataset. Using the exponential thinning method, the estimated volumes range from 0.91 to 17.66 times the full-dataset volume (1.15 to 22.7 km³), while the Weibull method results range from 0.75 to 2.89 times (0.96 to 3.56 km³). The largest estimates were observed only in the 10A dataset ( n i = 10 ), but they were not considered for further statistical treatment of data as they represent unrealistic estimates and are associated with volume estimate errors ≥250% (Fig. 5).

The exponential and Weibull methods show volume overestimations of ≥30% for smaller datasets. In the case of the exponential method, overestimation is mostly in datasets with ≤70 FCPs, while with the Weibull method overestimation is consistently observed with ≤50 FCPs. Conversely, underestimating the volume does not occur using exponential thinning. Using the Weibull method, however, does in some cases underestimate the volume relative to the full dataset (i.e., <20 %) typically when the dataset has n i = 60 - 100 .

For the exponential model, the p-value for the first and second hypotheses of interest are 0.3206 and 0.8867, respectively (Table 2). For the Weibull model, the p-values for the first and second hypothesis of interest are 0.4351 and 0.1923, respectively (Table 2).

Therefore, in both cases we have statistical support to guarantee that, for n i ≥ φ , both models provide similar results for the estimated tephra fall volume and the complete data. Specifically, the estimate for φ in the exponential case is 88.25 (95% confidence interval 27.54 - 148.96), and for the Weibull case it is 60.00 (95% confidence interval 38.89 – 81.11). Therefore, the preferred estimation for a sample size that guarantees a confident tephra fall volume is 89 and 60, with several isopachs ranging from 7 to 8 for the exponential and Weibull models, respectively (Table 2). However, if we consider a more conservative point of view (based on the upper limit of the 95% confidence intervals), the estimated sample sizes are 149 and 82 (resulting in 8-10 isopachs) for the exponential and Weibull models, respectively.

For datasets containing more than 90 FCPs, using the exponential model we find that the proximal isopachs return average volumes that contribute nearly 36% to the total volume, medial isopachs represent 52% and the distal isopachs are only 12% of the total volume.

Because of the range of volume estimates from the different datasets, we derive a range of eruption source parameters that differ up to one order of magnitude between the maximum and minimum values. For example, the total erupted mass (6.6 – 128.0 x10¹¹ kg) and peak MER (6.8 – 130.0 x 10⁶ kg s^-1) by the exponential method (Table 3), and similarly the erupted mass (5.3 – 20.3 x10¹¹ kg) and peak MER (5.8 – 21 x10⁶ kg s^-1) by the Weibull method (Fig. 6a). The minimum and maximum magnitude and intensity values vary by 1 unit for both volume estimate methods. Accordingly, the estimated column heights for these values are in the range 13.5-27.5 for the exponential method, and 12.8-17.7 for the Weibull method (Fig. 6b). Our best estimates for peak MER, M, I, and HT are 6.8±0.3 x 10⁶ kg s^-1, 4.82±0.02, 9.85±0.05, and 13.45±0.15 km, respectively, by considering the whole dataset. Comparatively, these values are similar to the minimum values estimated for ESPs reported in Table 3.

The approach chosen in our contribution aims to minimize bias in constructing isopach maps by limiting the use of knowledge on the eruption and its deposit, by using only the available data to construct the isopachs by hand-drawing knowing the position of the eruptive vent and the approximate dispersal direction (e.g., towards the E). However, this technique has its own flaws. For instance, Yang and Bursik (2013) identified a lack of adherence to the data in drawing isopachs for hand-drawn maps. Another limitation of hand-drawing is the influence of the local knowledge and experience of the geologist drawing the isopach contours, and the assumptions that modify isopach shapes (Engwell et al., 2015). These works have also proposed alternative strategies, such as the use of thickness trends and their separate modeling using segmented linear regression and ordinary kriging (Yang and Bursik, 2013) or the application of spatial statistical methods (Engwell et al., 2015). However, in recent ‘fully statistical’ algorithms for tephra volume estimation, expert judgment is unavoidable in making some choices (e.g., wind-direction, the limit of the tephra fall area, etc.; Rougier et al., 2022). Considering these differences and limitations, we can only compare the resulting volume estimates of both approaches as reported in the literature. For instance, the cubic B-spline curves under tension used by Engwell et al. (2015) yield volumes 20-40% smaller compared to hand-drawing isopachs, whereas the method used by Rougier et al. (2022) returned volumes that represent 37- 40% of the original bulk tephra estimate by Takarada and Hoshizumi (2020).

Bonadonna et al. (2015a) observed that varying the isopach contouring strategy might result in uncertainty up to 40% in determining associated erupted volume/mass. In comparison, the different volume integration strategies could result in uncertainties up to 70%, depending on the deposit exposure, data distribution, and eruption magnitude. From our dataset, we observe uncertainties up to 84% by comparing the volumes obtained from the exponential thinning and Weibull methods for datasets <60 FCPs, while the uncertainty is ≤41% for datasets ≥60 FCPs. In addition to these sources of uncertainty, Yang and Jenkins (2023) demonstrate that the misfit between the fitted curves and the isopach data as a result of model uncertainty could affect the tephra volume estimation; they propose that if half or more of the predictions from a certain fitted curve are outside the ±40% range, the curve should not be used for the volume calculation. We do not evaluate model vs. isopach misfit in this contribution.

The volume of the Puyehue-Cordón Caulle 2011-2012 tephra deposit has been studied based on proximal-medial data (0.5-1.0 km³, e.g., Amigo et al., 2012; Silva et al., 2012; Naranjo et al., 2017), and proximal-to-distal datasets (0.76-1.46 km³, e.g., Gaitán et al., 2011; Alloway et al., 2015; Pistolesi et al., 2015). The larger volumes estimated from proximal-to-distal data are consistent with the incorporation of larger datasets that include the outer limits of the deposit, and these results agree with our highest-confidence results (1.2-1.3 km³). This fact is also supported by the datasets and volume estimates in other recent Andean explosive eruptions. For instance, Watt et al. (2009) estimated a volume of 0.17 km³ for the Chaitén 2008-2009 eruption, mainly based on nine distal isopachs, while Alfano et al. (2011) also incorporated proximal-to-medial data and estimated 0.5-1.0 km³ utilizing a dozen isopachs. These results demonstrate the importance of considering medial (and proximal) data to constrain tephra volumes more accurately. The Calbuco 2015 tephra fall deposit provides another example. It was mapped in 57 proximal and medial locations by Castruccio et al. (2016), obtaining a range 0.28-0.38 km³, whereas the proximal-to-distal field observations of Romero et al. (2016) (65 sites) and Van Eaton et al. (2015) (163 sites) provided estimates of 0.27 km³ and 0.58±0.28km³, respectively. Later, Romero et al. (2021) used the full dataset of 150 field points to estimate a proximal-to-distal bulk volume in the 0.25-0.35 km³ range. All of them applied the Weibull, exponential and Power Law methods. The overall similarity of these estimates, in addition to a good distribution of the field data, agrees with the approximate number of field points required to achieve a 95% confidence interval using the Weibull method, as in this contribution. A remarkable example of high-resolution tephra mapping is observed for the Plinian eruption of Hudson volcano in 1991; Banks and Iven (1991) mapped a dozen proximal-to-distal isopachs and estimated a bulk volume of 4 km³, which does not differ considerably from Scasso et al. (1994) 4.35 km³ estimate, using fifteen isopachs and 200 field points.

Our results support findings from previous incompletely sampled tephra deposits that similarly found different sub-samples of the field data translates into uncertainties in erupted volume estimates, irrespective of eruption size (e.g., Andronico et al., 2014; Engwell et al., 2015; Buckland et al., 2020). Our exercise of the Puyehue-Cordón-Caulle 2011-2012 tephra fall deposit is a good example on the control of field points on accurately estimating tephra volumes and confirms previous evidence from Engwell et al. (2015), indicating that maps with <60 points produce volume errors generally >50%. Accordingly, small variations in the volume estimate are found at datasets >80-100 points. It is worth mentioning that Engwell et al. (2015) used a range of eruptions with published volumes from 0.003 to 8.2 km³ and a different isopach drawing method (cubic B-spline curves under tension). Also, because of the extensive thickness data available in our contribution, we have shown that integrating medial FCPs contributes to the largest portion of the total deposit volume (>50%), thus incompletely sampling deposits at medial areas will produce a larger uncertainty.

The careful investigation of tephra fall deposits with large field datasets provides the best means of determining accurate estimates of tephra volume, mass eruption rates, and interpretations of eruption dynamics. We determined that varying the number ( n i ) of FCPs, and consequently the supporting data for isopach drawing, may translate into order of magnitude ranges in the total mass and MER for 2011-2012 PCCVC eruption. Thus, we suggest authors be cautious when determining these parameters with limited FCPs. However, using field data collected at fewer locations can still be informative (Green et al. 2016), but has to consider huge uncertainties, sometimes greater than the volume estimate itself.

In light of our results and accordingly to the available literature, we suggest that tephra fall field mapping considers the total extent of the deposit: 1) if possible, covering from proximal to distal areas, but focusing on the proper characterization of proximal and medial areas as they account most of the erupted volume; 2) distributing FCPs isotropically, avoiding concentrating points along the crosswind or downwind axis only, which will probably result into a forced ellipsoidal distribution; 3) collecting, if possible, at least 60 to 90 field control points with these conditions, allowing the integration of 7-8 isopachs, and 4) estimating erupted volumes as a range, instead of absolute values (e.g., Constantinescu et al. 2022b). In this respect, we consider the best result is represented by the volume estimate derived from the integration of the largest number of FCPs. However, this choice is open to discussion and may depend on a series of physical and statistical considerations that we are not deciphering in detail. Considering the calculation of an exact volume is an unrealistic task, volcanologists are invited to provide a range of results considering different methods, and to report and discuss the different uncertainties that modify these results (e.g., Engwell et al. 2013; Bonadonna et al. 2015a; Biass et al. 2019). This is important as absolute estimates may represent a minimum or maximum, leading to underestimation or overestimation of ESPs, hazard mapping, or numerical modeling (e.g., Constantinescu et al. 2022b).