Introduction

Andesitic stratovolcanoes commonly produce complex explosive eruptions comprising multiple phases, styles, and intensities, from vent-clearing bursts and pyroclastic density currents (PDCs) to sustained sub-Plinian and Plinian columns (e.g., Cioni et al. 2003; Arce et al. 2005, 2012a; Carazzo et al. 2025; Suhendro et al. 2026). These transitions are driven by changes in magma chemistry, rheology, gas supply, magma mixing, conduit vesiculation and crystallization, and decompression and ascent rates, and are recorded in deposit architecture, componentry, juvenile textures, and geochemistry (e.g., Shea et al. 2010, 2012; Cassidy et al. 2018; Payet-Clerc et al. 2026). Reconstructing such transitions is essential for understanding eruption dynamics and defining more realistic hazard scenarios, especially at volcanoes whose largest eruptions have not occurred in historical time.

Pico de Orizaba (Citlaltépetl) is an active stratovolcano in eastern Mexico (Que-Salinas et al. 2025), the tallest in North America (5636 masl), with a long record of Pleistocene-Holocene andesitic-dacitic volcanism (Höskuldsson and Robin 1993; Carrasco-Núñez 2000). It last erupted in 1687 CE (Cantagrel et al. 1984), and about half a million people inhabit its flanks within 30 km of the summit crater (INEGI 2020), making it among Mexico’s most hazardous volcanoes. Recent work has shown that Holocene explosive activity at Pico was more frequent, stratigraphically complex, and process-diverse than previously recognized, involving repeated multi-phase eruptions with fall deposits and different types of PDCs (Sourisseau et al. 2025; Vásquez-Montoya et al. 2025; Arce et al. 2026). Among these, deposits formerly grouped as the “Citlaltépetl Pumice” (Rossotti et al. 2006) were recently redefined as two distinct early Holocene eruptions, the Lower and Upper Citlaltépetl, separated by different ages, paleosol development, contrasting juvenile assemblages, and different associated PDC units (Vásquez-Montoya et al. 2025). These eruptions provide a natural comparative pair to examine how two compositionally similar andesitic eruptions from the same volcano produced different depositional architectures, eruptive sequences, and hazard footprints; evaluate the implications of their contrasting behavior for explosive volcanism; and explore how differences in magma interaction and rheology may control transitions between fallout-dominated and PDC-dominated phases. Consequently, this work does not seek to re-establish the stratigraphic split, but takes the revised Lower-Upper framework as its starting point and focuses on the lithofacies architecture across multiple sectors of Pico de Orizaba, grain-size characteristics, componentry, petrography, microtextures, geochemistry, and quantification of source parameters of both eruptions separately, aiming to (1) refine and reconstruct eruptive sequences; (2) characterize juvenile components and magma signatures; and (3) unveil eruption styles and eruptive dynamics, and assess hazard implications. More broadly, this work contributes to understanding how Holocene andesitic stratovolcanoes can generate contrasting sub-Plinian to Plinian eruptions from the same edifice, and why this variability matters for long-term hazard assessment.

Geological background

The 650 ka-to-present Pico de Orizaba (Carrasco-Núñez 2000) marks the locus of active volcanism in the easternmost sector of the Trans-Mexican Volcanic Belt (TMVB) (Fig. 1a and b). The TMVB is a continental volcanic arc generated by subduction of the Cocos and Rivera plates beneath the North American plate along the Middle American Trench (Ferrari et al. 2012). In its eastern sector, volcanism is concentrated in a northeast-southwest- to north-south-trending range that includes the composite volcanoes Cofre de Perote, La Gloria, Las Cumbres, and Pico de Orizaba, plus the Sierra Negra satellite cone. The surrounding region also hosts diverse Quaternary small-volume volcanoes and associated lava fields (Carrasco-Núñez et al. 2021), all emplaced on a structurally complex basement composed mainly of folded and faulted northwest-southeast-trending Cretaceous limestones and shale successions of the Mesozoic Sierra Madre Oriental uplift (Fitz-Díaz et al. 2018; Díaz-Bravo et al. 2022).

Fig. 1
Fig. 1
Full size image

a Tectonic setting of the Trans-Mexican Volcanic Belt (light- to dark-gray areas), modified from Ferrari et al. (2012). Green lines represent depths of the subducting Rivera (RP) and Cocos plates. OGFZ: O’Gorman Fault Zone. Cities: Guadalajara (Gd), Toluca (To), Mexico City (Mx), Puebla (Pu), Xalapa (Xa), and Veracruz (Ve). Volcanoes: Mascota Volcanic Field (MVF), Colima Volcanic Complex (CVC), Michoacán-Guanajuato Volcanic Field (MGVF), Nevado de Toluca (NT), Iztaccíhuatl (IZ), Popocatépetl (P), Malinche (M), Cofre de Perote (CP), and Pico de Orizaba (PO). b Zoomed-in digital elevation model of the Pico de Orizaba (Citlaltépetl) edifice, summit crater at 5636 masl, and adjacent ring plain. The 85 study locations of this work are indicated and correspond to the exact same sites examined by Vásquez-Montoya et al. (2025). Location numbers shown in red font with bright-yellow halos correspond to Figs. 2 and 3. Closely clustered locations unresolved at map scale are grouped (e.g., 1–3). See the figure legend for further details

At edifice scale, Pico de Orizaba records repeated alternations between basaltic-andesitic to dacitic effusive and explosive activity, and between edifice growth and destruction through time (Höskuldsson and Robin 1993; Carrasco-Núñez 2000; Juárez-Arriaga et al. 2025), punctuated by two edifice-collapse debris avalanches at 290–210 ka (Carrasco-Núñez and Rose 1995) and ~16 ka (Carrasco-Núñez et al. 1993), and by formation of the modern Citlaltépetl cone and summit crater (Fig. 1b). Previous studies documented Holocene fallouts, block-and-ash and ash-and-scoria or pumice flows around the volcano, separated by repose intervals of 500–1500 years (Cantagrel et al. 1984; Höskuldsson and Robin 1993; Siebe et al. 1993; Carrasco-Núñez and Rose 1995; Rossotti et al. 2006). The Pleistocene-Holocene pyroclastic record has since been refined. Vásquez-Montoya et al. (2025) documented 40 distinct pyroclastic bedsets, each comprising interbedded fall and PDC deposits and separated by eruptive pauses of 100–800 years. At least 10 Pleistocene and 12 Holocene bedsets represent multi-phase explosive eruptions including onset dome-collapse flows, climactic eruption columns, and waning column-collapse PDCs. New work on the southern flank by Sourisseau et al. (2025) showed that at least ten pyroclastic flow units were emplaced between 8980 and 8170 BP, representing four Vulcanian eruptions and two dome-destruction events. Arce et al. (2026) described the ~5 ka Jamapa pyroclastic flow deposits on the northern flank, recording recharge-driven explosive behavior linked to ascent of andesitic magma from a deep (~20 km) reservoir into a shallower (~8 km) dacitic reservoir, where mingling/mixing preceded eruption.

The Lower and Upper Citlaltépetl

The < 8620 BP (9600–9500 cal BP) Lower Citlaltépetl and 8260 BP (9200 cal BP) Upper Citlaltépetl (14C dated and calibrated by Vásquez-Montoya et al. 2025) are two distinct explosive episodes formerly grouped within the legacy concept of the “Citlaltépetl Pumice” (Rossotti and Carrasco-Núñez 2004; Rossotti et al. 2006). In that framework, deposits now known to belong to separate eruptions were interpreted as parts of a single pumice-forming episode. However, recent stratigraphic revision showed that the succession actually records two discrete fallout-dominated eruptions, separated by intervening deposits and paleosol development, with contrasting 14C ages, juvenile clast characteristics (dominantly scoriaceous rather than pumiceous), and associated PDC units. According to Vásquez-Montoya et al. (2025), both Lower and Upper Citlaltépetl comprise a series of fallout and PDC layers produced during multiple eruptive phases (Fig. 2a). The Lower Citlaltépetl includes six layers, from onset PDCs (layers 1a–1b) to partially collapsing eruption columns (layers 2–4) and final column collapse (layer 5). The Upper Citlaltépetl comprises seven layers recording dilute PDCs (layers 1, 3, 5) interspersed with thick fallout layers (layers 2, 4, 6) derived from unsteady to steady eruption columns and capped by a massive column-collapse PDC (layer 7).

Fig. 2
Fig. 2
Full size image

a Composite column of deposits produced by the 9600–9500 cal BP Lower Citlaltépetl (LC) and 9200 cal BP Upper Citlaltépetl (UC) eruptions, modified from Vásquez-Montoya et al. (2025). Numbers indicate individual layers of each eruption. Red levels represent paleosols and/or reworked deposits separating the Lower and Upper eruptions. Grain-size classes are fine ash (FA), coarse ash (CA), fine lapilli (FL), coarse lapilli (CL), and blocks (B). b Volume-percentage (vol.%) componentry estimates for each layer, together with photographs of representative juvenile and non-juvenile clast types. c Grain-size distribution diagrams for all LC and UC layers except LC layer 1b. The first six diagrams correspond to the thickest and most distinctive fall deposits of UC layer 4 (UC4 in red contoured bars) and LC layers 2 and 4 combined (LC2+4 in orange bars), sampled at different locations (Loc. 04–56) within comparable distance from the summit crater. Inman (1952) mean (M) and sorting coefficient (σ) are indicated and correspond to UC4 and LC2+4 from left to right. The last three diagrams correspond to LC pyroclastic density current (PDC) deposits (layers 1a, 3, and 5), UC PDC deposits (layers 1, 3, 5, and 7), and secondary UC fallout layers 6 and 2. d Walker (1971) mean-versus-sorting diagram used to classify deposits as fallouts or PDCs

Methods

Fieldwork, stratigraphy, and mapping

Eighty-five locations exposing Lower and Upper Citlaltépetl deposits were examined 3.5–25 km from the summit crater (Fig. 1b, Online Resource 1). Detailed lithostratigraphic descriptions were completed at every location and included all layers recognized by Vásquez-Montoya et al. (2025) (Fig. 2a). Here, “layer” refers to an individual depositional unit emplaced during a single fallout or PDC event, bounded by locally erosive, gradational, or sharp contacts produced by changes in grain size and/or framework. The broader stratigraphic architecture and bounding relationships of the Lower and Upper Citlaltépetl follow the framework of Vásquez-Montoya et al. (2025). At each location, we recorded layer contacts, thickness, geometry, grain size (White and Houghton 2006), framework, grading, juvenile and lithic component characteristics (Fig. 2b), and syn-depositional deformation structures, including impact sags, and diffuse and cross-bedding where present (e.g., Arce et al. 2012a; Cas et al. 2008, 2023). These observations were used to identify and correlate fallout and PDC layers, construct stratigraphic logs, define lithofacies from recurrent physical attributes and internal structures (lithofacies descriptions in Table 2, Appendix), and recognize lithofacies transitions recording proximal-to-distal depositional changes and lateral variations perpendicular to the main dispersal axes (e.g., Torres-Orozco et al. 2018, 2025).

Isopach and isopleth maps were drawn digitally in ArcMap 10.5 from deposit thicknesses and maximum clast sizes measured from the three orthogonal axes of the five largest juvenile scoria/pumice and five largest dense clasts at each location. Fallout isolines were interpolated using combined kriging and natural neighbors (Torres-Orozco et al. 2017a; 2025) (Online Resource 1). PDC deposits were interpolated using kriging and inverse-distance weighting because PDC data points were evenly distributed and less strongly clustered across the study area (cf. Plank et al. 2023). Incomplete fallout contours were geologically corrected using adjacent better-constrained isolines and extrapolated dispersal trends following elliptical shape factors ≤ 1, aspect ratios ≥ 1, and eccentricities of 0–1.05 (e.g., Cioni et al. 2000; Sulpizio 2005; Klawonn et al. 2014).

Granulometry and componentry

Thirty-five bulk samples (1.5–3 kg each) were collected from all Lower and Upper Citlaltépetl layers. Samples were dry-sieved at 1 ϕ intervals (from −5 to > 4 ϕ). Grain-size distributions and Inman (1952) statistical parameters were estimated (Fig. 2c) (Online Resource 1). Mean grain size (Mϕ) and sorting coefficient (σ1ϕ) were used to characterize the deposits following Walker (1971) (Fig. 2d). All samples were further divided into representative coarse (−2 to −3 ϕ) and fine-grained (0 to 1 ϕ) fractions, each consisting of 100 clasts; thus, 400 clasts per sample were counted manually under a binocular stereomicroscope. Clasts were classified after White and Houghton (2006) and Cas et al. (2023) into juvenile (scoria, banded scoria, pumice, dense lava) and non-juvenile types (red-oxidized accessory clasts and some heavily altered whitish lithics) (Fig. 2b). Dense lava clasts were further subdivided into lava type 1 (gray porphyritic lava) and lava type 2 (dark gray-black, porphyritic vitric lava). Component abundances were estimated for all layers except Lower Citlaltépetl layer 1b, which is rarely preserved (complete data in Online Resource 1).

Petrographic and chemical analyses

Petrographic analysis of 20 juvenile clast thin sections (11 Lower and nine Upper Citlaltépetl) used plane- and cross-polarized light microscopy to characterize mineral assemblages, crystal textures and sizes, vesicularity, crystallinity, and domain relationships. Vesicularity and crystallinity were estimated by manually segmenting and counting vesicles and crystals within 6.25 mm2 per thin section. Treating each analyzed area as 100 vol.%, vesicle and crystal area fractions were converted to vol.%. All samples were analyzed using comparable areas and criteria; however, replicate measurements were not performed. Following Zellmer (2021), microlites were < 30 µm wide, and phenocrysts comprised micro-/mesocrysts (30–300 µm long), macrocrysts (500–1000 µm), and megacrysts (> 1000 µm). Vesicles were classified as spherical, subspherical, oblate, coalesced, tubular, or fibrous (sheared) (Torres-Orozco et al. 2023, 2024). Juvenile clasts were classified as dense-incipiently vesicular (0–20 vol.% vesicles), poorly vesicular (20–40 vol.%), vesicular (40–60 vol.%), and highly vesicular (> 60 vol.%) (Houghton and Wilson 1989). Microvesicles were < 150 µm in diameter, and high microvesicularity indicated that they comprised ≥ 50 vol.% of total vesicularity normalized to 100 vol.% (Klug et al. 2002; Shea et al. 2010; Torres-Orozco et al. 2025). Major- and trace-element compositions were measured from 28 individually selected and powdered juvenile scoria, banded scoria, pumice, and dense lava clasts using a Rigaku ZSX Primus II sequential X-ray fluorescence spectrometer at the Laboratorio Nacional de Geoquímica y Mineralogía, Universidad Nacional Autónoma de México. These comprised 11 Lower and 10 Upper Citlaltépetl analyses and seven reference analyses from other Pico de Orizaba deposits. Complete dataset is available in Online Resource 1.

Eruption source parameters

Fall deposit volumes were estimated from isopach data using exponential (Fierstein and Nathenson 1992; Pyle 1989, 2015), proximal extrapolation (Sulpizio 2005), and Weibull-fit methods (Bonadonna and Costa 2012, 2013). Most estimates were consistent within their corresponding standard deviations. PDC deposit volumes were calculated from interpolated thickness-interval polygons constructed in ArcMap (Online Resource 1). Overlapping areas were subtracted to obtain the net area of each interval, and volumes were calculated using the mean thickness of each range and summed to obtain total PDC volumes (Arce et al. 2005; Torres-Orozco et al. 2017a). Dense-rock-equivalent (DRE) volumes were calculated by multiplying the fallout and PDC bulk volumes by the ratio of representative bulk deposit densities to an assumed dense-rock andesite magma density of 2.5 g cm−3 (e.g., Al-Shammary et al. 2018; Torres-Orozco et al. 2025). Eruption column heights were estimated from isopleth data following Carey and Sparks (1986), Sulpizio (2005), and Bonadonna and Costa (2012, 2013) (Online Resource 1). Volume and mass eruption rates (MERs) were calculated following Sparks (1986), Carey and Bursik (2000), and Mastin et al. (2009). Eruption mass and magnitude were estimated following Wilson (1976) and Pyle (2015). Eruption styles were classified using deposit-thinning rates of Pyle (1989) and Bonadonna and Costa (2013), and the classification scheme of Mastin et al. (2009). Complete data are shown in Table 1.

Table 1 Estimated fallout-layer and total-deposit volumes, and eruption source parameters
Full size table

Results

Stratigraphy, lithofacies, and grain size

Lithofacies correlations show that both Lower and Upper Citlaltépetl are best preserved in the eastern and southeastern sectors, where the thickest and most complete successions occur, whereas northward and south-southwestward they thin, become less continuous, and are commonly affected by reworking or erosion (Figs. 3 and 4). Six facies groupings are recognized: distal-medial northern, medial-proximal northeastern, proximal eastern or southeastern, and either medial-proximal or medial-distal southern-southwestern (Fig. 4). Across these sectors, fallout layers are the most laterally persistent and are generally represented by clast-supported lapilli lithofacies, bearing scoria or scoria and pumice clasts (lithofacies descriptions in Table 2, Appendix). In contrast, PDC layers are dominated by matrix-supported tuff lithofacies (massive to diffusely stratified or cross-bedded), and composed of scoria, scoria and pumice, or lava and lithics (Figs. 3 and 4). The Lower Citlaltépetl is less well preserved than the Upper Citlaltépetl and commonly survives as discontinuous remnants in northern and southwestern sectors (Figs. 3f and g and 4a and f). Its basal layers 1a–1b are represented by massive tuffs. Fallout layers 2 and 4 form the thickest part of the sequence and are dominated by clast-supported scoria-pumice lapilli, whereas layers 3 and 5 are commonly expressed by massive tuffs (Figs. 3 and 4). The Upper Citlaltépetl is more continuous and widespread and is dominated by a thick marker layer 4 composed mainly of clast-supported scoria lapilli. Layer 2 is distinctive for its high dense lava and lithic contents, whereas layers 5–7 show the greatest lithofacies variability, including alternations of massive and cross-bedded tuffs and massive lapilli (Fig. 3a–e). Together, these correlations show that lithofacies vary not only between sectors, but also within individual transects (Figs. 1b and 4).

Fig. 3
Fig. 3
Full size image

Photographs of Lower (LC) and Upper Citlaltépetl (UC) deposits exposed at different locations. a–b represent proximal deposit lithofacies, c–e correspond to medial lithofacies, and f–g represent medial-distal and/or lateral lithofacies variations. Lithofacies are indicated by alphabetic codes inside white boxes (full lithofacies nomenclature in Table 2, Appendix), whereas individual layers are indicated by numbers. Yellow lines mark erosive or sharp contacts defined by paleosols, weathered ash (WA), or reworked deposits (Rw), separating the Lower and Upper eruptions. Solid and dashed black lines represent sharp or locally erosive, and gradational contacts, respectively, marked by changes in internal deposit framework and/or grain size. Radiocarbon ages in (c) and (e) were determined by Vásquez-Montoya et al. (2025): (*1) 9221 cal BP, (*2) 9229 cal BP, (*3) 9605 cal BP

Fig. 4
Fig. 4Fig. 4
Full size image

Stratigraphic-log correlations and lithofacies transitions of the Lower (LC) and Upper Citlaltépetl (UC) deposits along six transects (a–f), arranged clockwise from the northern to southern flanks of the volcano, as indicated in Fig. 1b. For complete lithofacies nomenclature, see Table 2 (Appendix). Interpretation of the stratigraphic logs, including major LC and UC contacts (red levels), clast and deposit colors and framework, grain size, and layer contacts, is provided in the figure’s legend. Overall, deposits increase in thickness from northern to eastern and southeastern locations and then decrease again toward the south and southwest. Accordingly, six lithofacies groups were distinguished: a distal-medial northern facies, b medial-proximal northeastern facies, c–e proximal downwind eastern to southeastern facies, and (f) proximal to medial-distal southern-southwestern facies. Facies were further associated with downwind and crosswind directions based on the isopach and isopleth maps. In addition, deposits within individual transects and facies groups also change with increasing distance from the summit crater of Pico de Orizaba: fallout layers tend to remain relatively consistent along transects, whereas pyroclastic density current layers are variable, commonly pinching-and-swelling or even disappearing between locations (e.g., LC layer 5 from locations 64–65 to 61 in transect d)

The grain-size distributions of the thickest and most representative fallout layers, Upper Citlaltépetl layer 4 (UC4) and Lower Citlaltépetl layers 2 and 4 (LC2+4), were analyzed at six locations spanning northern to southern-southwestern sectors clockwise, where they correspond to different lithofacies variants (Figs. 2c and 3). Despite these spatial changes, their grain-size distributions remain broadly consistent. All UC4 and LC2+4 samples, and Upper layers 2 and 6, are unimodal and plot in the fallout field of the Walker (1971) mean-versus-sorting diagram (Fig. 2d), whereas Lower layers 1a, 3, and 5, and Upper layers 1, 3, 5, and 7, range from multimodal to unimodal and negatively skewed distributions, and plot in the PDC field. These grain-size data support the layer-by-layer stratigraphic framework of Vásquez-Montoya et al. (2025) (Fig. 2a).

Componentry

The Lower Citlaltépetl is characterized by a heterogeneous juvenile assemblage (Fig. 2b). Its thickest layers 2–4 are dominated by banded scoria (37–52 vol.%). Basal layer 1a is exceptional in comprising 88 vol.% dense lava clasts (types 1 and 2), whereas layer 5 contains nearly balanced proportions of scoria (35 vol.%) and pumice (21 vol.%), together with subordinate banded scoria (17 vol.%), dense lava (20 vol.%), and lithics (7 vol.%). The Upper Citlaltépetl layers 3–7 are predominantly scoriaceous (52–79 vol.%), with pumice absent and banded scoria present only in minor amounts in layers 3–4 (4–7 vol.%). Basal layers 1–2 stand out for containing 59–69 vol.% lava type 1 and the highest lithic contents (12–16 vol.%).

Petrography and microtextures

In thin section, Lower Citlaltépetl banded clasts show direct juxtaposition of scoria and subordinate pumice domains with contrasting vesicularity and mineral assemblages (Fig. 5a–d). Scoria is highly vesicular (~55–62 vol.%) and microvesicular (> 50 vol.% of total vesicularity), containing abundant spherical to oblate microvesicles (< 5–10 µm in diameter) and fewer larger oblate vesicles (~100–250 µm) (Fig. 5e). Its mineral assemblage (Fig. 5e–i) is dominated by macro- to megacrysts of euhedral amphibole (40 vol.%) and clinopyroxene (25 vol.%), followed by subhedral to euhedral plagioclase meso- to megacrysts (30 vol.%), and minor Fe-Ti oxide microcrysts (5 vol.%). Plagioclase commonly shows sieve textures and zoning. In contrast, pumice is vesicular (~45 vol.%), and is characterized by large coalesced and oblate vesicles (~100–250 µm), together with fewer elongated and coalesced microvesicles (~50 µm in diameter) (Fig. 5j and k). Its assemblage (Fig. 5l and m) comprises microlites and meso- to megacrysts of euhedral to subhedral plagioclase (60 vol.%) and orthopyroxene (30 vol.%), with subordinate Fe-Ti oxides (10 vol.%). Plagioclase occurs as acicular groundmass microlites, small tabular crystals, and larger sieve-textured crystals.

Fig. 5
Fig. 5
Full size image

Thin-section photomicrographs in cross- and plane-polarized light showing minerals and groundmass microtextures of Lower Citlaltépetl juvenile clasts. a–d Banded scoria clasts showing juxtaposition of scoria and pumice domains and abundant amphibole crystals. e–i Highly vesicular and microvesicular scoria clasts and minerals. Individual vesicles in blue color (e′) were manually segmented and counted for modal and shape analyses. The same procedure was followed for modal crystal analysis. j–m Vesicular pumice clasts and minerals. Individual vesicles in blue color (k′) were manually segmented and counted for modal and shape analyses. Insets e′ and k′ only exemplify the procedure applied to total 6.25-mm2-thin section areas for modal vesicle and crystal analyses

The Upper Citlaltépetl is dominated by scoria, and most banded clasts are produced by textural heterogeneity within the same scoria, involving juxtaposition of microvesicular scoria 1 (Sc1), crystalline scoria 2 (Sc2), and sheared scoria 3 (Sc3), the latter resembling a deformed version of Sc1 (Fig. 6a–f). Juvenile pumice is absent, except for rare highly oxidized pumice-like domains preserved within some banded clasts (Fig. 6a and b). Sc1 is highly vesicular (~57–66 vol.%) and microvesicular (~60–70 vol.% of total vesicularity), containing mostly subspherical microvesicles (2–10 µm in diameter) and a second population of ~100–150 µm (Fig. 6c–f). Sc2 is vesicular (~39–45 vol.%), highly microcrystalline (~63–69 vol.% groundmass microlites and microcrysts), and characterized by mainly oblate vesicles (~100–500 µm) (Fig. 6g and h). Sc3 is highly vesicular (~62–68 vol.%) and dominated by tubular to fibrous vesicles < 10 µm in diameter but up to ~500 µm long, producing a sheared appearance (Fig. 6b, e, and f). The mineral assemblage (Fig. 6i–t) is dominated by subhedral clinopyroxene meso- to macrocrysts, commonly forming aggregates (24–34 vol.%), followed by euhedral orthopyroxene meso- to macrocrysts (12–19 vol.%), subhedral amphibole micro- to macrocrysts (5–15 vol.%), euhedral to anhedral plagioclase microlites and micro- to megacrysts (28–41 vol.%), subhedral Fe-Ti oxide microlites and microcrysts (7–12 vol.%), and rare subhedral olivine mesocrysts (2 vol.%). Plagioclase occurs as acicular microlites, tabular micro- to mesocrysts, zoned crystals with dissolution features, and sieve-textured megacrysts.

Fig. 6
Fig. 6
Full size image

Thin-section photomicrographs in cross- and plane-polarized light showing minerals and groundmass microtextures of Upper Citlaltépetl juvenile clasts. a–f Banded scoria clasts showing rare juxtaposition of scoria and highly altered (oxidized) pumice, and predominant juxtaposition of contrasting scoria types, ranging from highly vesicular, microvesicular scoria 1 (Sc1), through highly microcrystalline, poorly vesicular scoria 2 (Sc2), to highly vesicular, sheared scoria 3 (Sc3). Note differences in vesicle shapes: subspherical (Sc1), oblate (Sc2), and tubular to fibrous (Sc3). g and h Individual vesicles in blue color were manually segmented and counted for modal and shape analyses. The same procedure was employed for modal crystal analysis. The segmented blue vesicles in f and h only exemplify the procedure applied to total 6.25-mm2-thin section areas. i–t Mineral assemblage characteristic of all scoria clasts. Both clinopyroxene aggregates and individual crystals are common. Rare olivine crystals were also identified. Pervasive plagioclase macro- to megacrysts with sieve textures and zoning coexist with euhedral, acicular to tabular plagioclase microlites and micro- to mesocrysts

Geochemistry

The chemical compositions of Lower and Upper Citlaltépetl juvenile clasts range from basaltic-andesite to high-Si andesite and define a medium-K sub-alkaline trend (Fig. 7a). Most Upper Citlaltépetl samples cluster tightly within the andesitic field, between ~58 and 60 wt.% SiO2. In contrast, the Lower Citlaltépetl shows a broader compositional spread, from ~56 to 62 wt.% SiO2. Lower pumice defines the most evolved compositions of the dataset and is systematically depleted in Fe2O3, MgO, TiO2, MnO, CaO, and P2O5, while showing comparatively high Na2O and SiO2 values (Fig. 7a–c and Online Resource 1). By contrast, associated scoria, dense lava, and banded scoria partly overlap the Upper field, but collectively extend toward slightly more mafic and alkaline compositions, particularly in MgO, TiO2, CaO, and P2O5 (Fig. 7b and c and Online Resource 1). Trace-element ratios and concentrations also distinguish both eruptions (Fig. 7d–l). The Upper Citlaltépetl forms a relatively restricted compositional cluster, with limited dispersion in Sr/Y, Rb-Sr, Zr/Nb, Zr/Y, and Nb/Y. In contrast, the Lower Citlaltépetl is more dispersed. Its pumice has the lowest Y contents, but the highest Zr/Nb, Th/Nb, Sr/Y, Zr/Y, and Nb/Y values, whereas the remaining components define a more extended compositional range (Fig. 7g–i). The Th/Nb versus Zr/Nb plot is particularly diagnostic (Fig. 7e): except for Lower pumice, which extends toward higher values in both ratios, the Upper and Lower Citlaltépetl largely overlap in Th/Nb, whereas the Upper is clearly shifted toward higher Zr/Nb values. In Ni- and Cr-based plots (Fig. 7j and k), the Upper Citlaltépetl occupies a relatively restricted intermediate field, whereas Lower pumice defines the low-compatible-element end member and associated scoria, dense lava, and banded scoria extend toward higher Ni and Cr values (Fig. 7a–l and Online Resource 1).

Fig. 7
Fig. 7
Full size image

Whole-rock major- and trace-element compositions (anhydrous basis) of the Lower and Upper Citlaltépetl juvenile components. a Total alkalis versus silica diagram (Le Maitre et al. 2002), with the sub-alkaline/alkaline division of Irvine and Baragar (1971). Reference data sources are (1) Vásquez-Montoya et al. (2025), (2) Rossotti et al. (2006), (3) Carrasco-Núñez and Rose (1995), (4) Carrasco-Núñez et al. (2010), and (5) Rodríguez et al. (2002). b and c Major-element binary diagrams for MgO and TiO2. Remaining major-element binary diagrams are provided in Online Resource 1. d–i Trace-element concentration and ratio binary diagrams

Dispersal, eruptive volumes, and column height

Lower Citlaltépetl PDC layers show broad proximal-medial spreading toward the east, but also important northward and southward components (Fig. 8a–c). Layer 1a is especially notable for its strongest interpolation toward the south-southwest, where thickness reaches 77 cm (Fig. 8a). Layer 5 shows a comparable pattern, but with the greatest interpolated thickness (maximum 300 cm) concentrated in a restricted southeastern sector (Fig. 8c), whereas layer 3 is much smaller and radially distributed, with comparatively uniform interpolated thicknesses of < 6 cm (Fig. 8b). In contrast, the main fallout layers 2+4 were dispersed southeastward, as indicated by 5–20 cm isopachs covering 100–500 km2 (Fig. 8d). Corresponding isopleths also define a southeastern trend (Fig. 8e and f). Estimated bulk PDC volumes are 0.014 km3 (layer 1a), 0.006 km3 (layer 3), and 0.05 km3 (layer 5). Fallout volume is 0.13 ± 0.01 km3 (0.07 km3 DRE), giving a total Lower Citlaltépetl bulk volume of 0.20 km3 (0.10 km3 DRE; Table 1).

Fig. 8
Fig. 8
Full size image

Distribution maps of fall and pyroclastic density current (PDC) deposits of the Lower Citlaltépetl (LC). White-contoured numbers indicate thicknesses or average clast sizes measured in the field at locations where deposits are exposed. a–c Thickness maps of LC PDC layers 1a, 3, and 5, interpolated using combined kriging and inverse-distance weighting. The PDC thickness legend in panel a applies to panels a–c. d Isopach map of fallout layers 2+4. Yellow dashed squares delimit the mapped areas of PDC layers 3 and 5 for comparison with the total fallout extent. e and f Isopleth maps of dense clasts (lava and/or lithics) and scoria from fallout layers 2+4, interpolated using combined kriging and natural neighbors; the same method was used to define the isopach contours. The isopleth and clast-size legends in panel e apply to panels e and f. Complete data are provided in Online Resource 1

Upper Citlaltépetl PDC layers are the most widespread of the two eruptions. Layer 1 was distributed mainly northward and southward, reaching up to 16 cm in northern gullies (Fig. 9a), whereas layers 3 and 5 spread across the broadest eastern areas mapped for either eruption, although interpolated thicknesses remain modest, generally only up to 10–20 cm in the thickest restricted eastern sectors (Fig. 9b and c). Fall deposits were dispersed mainly eastward and northeastward (Fig. 9d–f). Layer 2 defines a northeastern dispersal axis, with 4–8 cm isopachs covering 200–650 km2, and < 0.5 cm scoria and lava/lithic clasts dispersed to 25 km downwind (Fig. 10a and b). Layer 4 is the largest, thickest, and most widely dispersed fallout unit of both eruptions, with 5–30 cm isopachs covering 160–1500 km2 eastward (Fig. 9e), and 1 cm clasts reaching 25 km downwind (Fig. 10c and d). Layer 6 is more restricted, with 5–20 cm isopachs covering 60–430 km2 and isopleths extending 16–19 km downwind (Figs. 9f and 10e and f). Bulk PDC volumes are 0.009 km3 (layer 1), 0.02 km3 (layer 3), and 0.05 km3 (layer 5). Average fallout volumes are 0.10 ± 0.02 km3 (layer 2), 0.32 ± 0.04 km3 (layer 4), and 0.09 ± 0.01 km3 (layer 6), yielding a total Upper Citlaltépetl bulk volume of 0.57 km3 (0.30 km3 DRE; Table 1). Corresponding average column heights are 21 ± 3.8 km (layer 2), 29 ± 0.9 km (layer 4), and 24 ± 0.4 km (layer 6), whereas the Lower Citlaltépetl average column height is 25 ± 0.9 km (layers 2+4) (Fig. 11 and Table 1).

Fig. 9
Fig. 9
Full size image

Distribution maps of fall and pyroclastic density current (PDC) deposits of the Upper Citlaltépetl (UC). White-contoured numbers indicate thicknesses measured in the field at locations where deposits are exposed. a–c Thickness maps of UC PDC layers 1, 3, and 5, interpolated using combined kriging and inverse-distance weighting. The PDC thickness legend in panel a applies to panels a–c. d–f Isopach maps of fallout layers 2, 4, and 6, interpolated using combined kriging and natural neighbors. For comparison with fallout extents, yellow dashed squares delimit the mapped areas of PDC layers 1, 3, and 5, which preceded the eruption columns that produced layers 2, 4, and 6, respectively. The isopach legend in panel d applies to panels d–f. Complete data are provided in Online Resource 1

Fig. 10
Fig. 10
Full size image

Isopleth maps of dense clasts (lava and/or lithics) and scoria from the Upper Citlaltépetl (UC) fallout layers 2 (a and b), 4 (c and d), and 6 (e and f), interpolated using combined kriging and natural neighbors. White-contoured numbers indicate average clast sizes measured in the field at locations where the deposits are exposed. Note that each panel has its own clast-size legend. Complete data are provided in Online Resource 1

Fig. 11
Fig. 11
Full size image

Lower (LC) and Upper Citlaltépetl (UC) eruption-style classifications and source parameters. a and b Eruption-style classification based on deposit-thinning rates and maximum dense-clast size, following a Walker (1973) and Pyle (1989, 2015), and b Bonadonna and Costa (2013). c Dense-rock-equivalent (DRE) eruptive volume versus column height for individual fallout layers and total LC and UC eruption deposits. Eruption magnitudes (M) from 4 to 6 were calculated following Pyle (2015) using the equation log10M [in kg] − 7. d Eruption-style classification based on volume eruption rate and mass eruption rate (MER) versus column height, modified from Mastin et al. (2009), using final average MER values from Mastin et al. (2009), Sparks (1986), and Carey and Bursik (2000). References: (1) Carey and Sigurdsson (1986), (2) Mastin et al. (2009), (3) Cas et al. (2023), (4) Saucedo et al. (2010), (5) Bonadonna and Costa (2012, 2013), (6) Torres-Orozco et al. (2017a, 2017b), (7) Carazzo et al. (2012), (8) Carey and Sigurdsson (1987), Macedonio et al. (1988), (9) Arce et al. (2005), (10) Avellán et al. (2014), (11) Sosa-Ceballos et al. (2012), (12) Capra et al. (2006), (13) Siebe et al. (2017), (14–15) Arce et al. (2012a, 2012b), (16) Torres-Orozco et al. (2025). Eruption source parameters are presented in Table 1. Complete data are available in Online Resource 1

Eruption style, intensity, and magnitude

Based on the thinning-rate parameters of Pyle (1989) and Bonadonna and Costa (2013), the Lower and Upper Citlaltépetl fallout layers correspond to sub-Plinian to Plinian eruptions (Fig. 11a and b). Within this range, Upper Citlaltépetl layer 4 occupies the highest-intensity end member, whereas layer 2 plots as the least intense, and layer 6 and Lower Citlaltépetl layers 2+4 remain intermediate. The same pattern is reproduced in DRE volume versus column-height space and in volume eruption rate and MER space (Fig. 11c and d), where layer 4 plots in the Plinian field and the remaining layers cluster between sub-Plinian and Plinian styles. Average mass eruption rates range from 3.3 × 107 to 3.5 × 108 kg s−1 for Upper layers 2 and 4, respectively, with intermediate values of 1.1 × 108 kg s−1 for layer 6 and Lower layers 2+4 (Table 1). Following Wilson (1976) and Pyle (2015), erupted masses for individual fallout layers range from 1.1 × 1011 to 4.2 × 1011 kg, and increase to 2.6 × 1011 and 7.5 × 1011 kg for the total Lower and Upper Citlaltépetl, respectively. These yield Upper Citlaltépetl fallout-layer magnitudes of 4.0–4.6, whereas inclusion of mapped PDC deposits raises total eruption magnitudes to 4.4 and 5.0 for the Lower and Upper Citlaltépetl, correspondingly (Fig. 11c).

Discussion

The former “Citlaltépetl Pumice” succession of Rossotti et al. (2006) was split into two eruptions by Vásquez-Montoya et al. (2025). The present work further supports that separation and shows that several key units were previously misinterpreted in eruptive style, correlation, and dispersal. In the Lower Citlaltépetl (LC), Rossotti’s fallout layers A0 and B–C are here redefined as PDC deposits (LC layers 1a and 5, respectively). Their Plinian layer A corresponds to our main sub-Plinian column phase (LC2+4), which was partially collapsing to generate LC layer 3. In the Upper Citlaltépetl (UC), several Rossotti layers (2, 3, 4) do not belong to the eruption, in agreement with Sourisseau et al. (2025) and Vásquez-Montoya et al. (2025), whereas the remaining units record a northeastward early sub-Plinian fallout (UC2), a broader PDC phase (UC3), a wider eastward-dispersed Plinian fallout (UC4), a highly energetic PDC (UC5), and two final sub-Plinian and PDC phases (UC6–7), none of which was properly resolved before. Further details and evidence are discussed below.

The 9600–9500 cal BP Lower Citlaltépetl sub-Plinian eruption

The Lower Citlaltépetl eruption began with destruction of a summit dome, before rapidly evolving into a sub-Plinian fallout-producing phase, punctuated by dilute PDCs and ending with a column-collapse event (Figs. 8 and 12a–d). Its total eruptive volume (0.1 km3 DRE) and magnitude (4.4) place it among the most explosive early Holocene events of Pico de Orizaba (Fig. 11). The initial stage is represented by layer 1a, a thin deposit dominated by lava- and lithic-rich ash-tuff lithofacies, with scarce pumice, diffuse to massive internal structures, and only slightly erosive basal contacts (Figs. 3, 4, and 12b). Its interpolated thickness map shows a strongly localized PDC axis toward the south-southwest (Fig. 8a), which is most consistent with emplacement by a block-and-ash flow (BAF) produced by initial summit-dome collapse from within the Pico de Orizaba crater (Fig. 12a and b) (Macorps et al. 2018; Capra et al. 2018, 2024). We infer that dense channel-confined BAF facies were largely removed by later erosion, and the preserved record corresponds mainly to associated dilute PDCs.

Fig. 12
Fig. 12
Full size image

Sketch of inferred eruptive phases, hazards, and first-order eruptive dynamics of the Lower and Upper Citlaltépetl eruptions. Lower Citlaltépetl: a possible pre-eruptive summit-dome/lava-flow effusion-growth stage. b layer 1a records initial summit-dome collapse, generating a block-and-ash-flow (BAF) avalanche and associated ash-cloud PDCs. Minor silicic-andesite extraction is inferred from the presence of subordinate pumice in the BAF deposits. c Waxing conduit/vent-clearing phase followed by an unsteady, partially collapsing sub-Plinian column that produced fallout layers 2+4 and associated PDC deposits (layers 1b and 3). Magma mixing reached its maximum, generating the highest banded scoria contents. d layer 5 records complete column collapse and concentrated massive PDCs. Magma mixing waned, and deposits became dominated by individual scoria and pumice clasts. Upper Citlaltépetl: e layer 1 records dense-clast-rich vent/conduit-clearing PDCs, followed by an early steady sub-Plinian column (layer 2) fed by mafic-andesite magma. f Layer 3 records a highly mobile PDC, followed by a sustained steady Plinian column that produced scoria fallout layer 4. Faint mingling may have generated minor banded scoria-pumice clasts, but the dominant control was a well-connected conduit and gas-rich mafic-andesite magma with strong internal rheological heterogeneity. g Progressive destabilization of the Plinian column produced broad PDC deposits (layer 5) and a waning, unsteady sub-Plinian column (layer 6). This phase likely reflects tapping of progressively outgassed magma, conduit widening, and waning gas-rich magma supply. h Layer 7 records final complete column collapse, producing concentrated massive PDCs

The eruption transitioned into a conduit-clearing stage represented by layer 1b, a poorly preserved tuff layer that shares the same general PDC character as layer 1a, but differs markedly in clast composition. Layer 1b contains mainly banded scoria and pumice of ash to fine-lapilli size, with much less dense lava and lithic clasts, indicating that juvenile magma was already being erupted (Figs. 3, 4, and 12c). Because layer 1b is immediately followed by the main fall deposit (layers 2+4) and is compositionally similar to it, it most likely records waxing unstable eruptive conditions just before establishment of the main column (Branney and Kokelaar 2002; Kim et al. 2014). The climactic phase generated the principal fall deposits of the eruption (layers 2+4), composed mainly of clast-supported lapilli lithofacies and dominated by banded scoria (Figs. 2b, 3, and 4). Their grain-size distributions are remarkably consistent across contrasting lithofacies sectors, and all samples plot in Walker’s fallout field (Fig. 2c and d), supporting interpretation as products of a single main fallout. This fallout was dispersed southeastward (500 km2) (Fig. 8d–f). Column-height estimates of 25 ± 0.9 km, together with average MER values of 1.1 × 108 kg s−1, indicate a sustained sub-Plinian column (Figs. 11 and 12c and Table 1). The thin layer 3, intercalated between layers 2 and 4, is represented by diffusely bedded tuff lithofacies, and its small, restricted, nearly radial distribution, with interpolated thicknesses of < 10 cm (Fig. 8b), indicates a comparatively minor dilute PDC deposited during temporary column unsteadiness (Cioni et al. 2003, 2008; Torres-Orozco et al. 2018).

The final stage is represented by PDC layer 5, dominated by massive tuff lithofacies and typically showing sharper basal contacts than the preceding PDC layers, locally eroding layers 2–4 (Figs. 3 and 4). Its clast population remains dominated by scoria, pumice, and banded scoria, with no strong increase in dense lava or lithic content, suggesting that the change from fallout to PDC deposition did not result from progressive densification of erupted material, but possibly from increased MER (Carazzo et al. 2025). Layer 5’s highest volume (0.05 km3), together with its mapped distribution and interpolation, dispersed mainly southeastward, where it reached maximum thicknesses in restricted locations (Figs. 4 and 8c), indicates a more concentrated and confined PDC, interpreted as the product of total eruption-column collapse (Branney and Kokelaar 2002; Shea et al. 2012; Scarpati et al. 2020), marking the end of the Lower Citlaltépetl episode (Fig. 12d).

The 9200 cal BP Upper Citlaltépetl sub-Plinian to Plinian eruption

The Upper Citlaltépetl began with vent-clearing activity, followed by successive sub-Plinian to Plinian fall phases, and ending with progressive destabilization, partial collapse, waning fallout, and final complete column collapse (Figs. 9, 10, 11, and 12e–h). Its total eruptive volume (0.3 km3 DRE) and magnitude (5.0) make it the largest and most intense event documented in this work (Fig. 11 and Table 1). The initial stage is represented by layer 1, which closely resembles Lower layer 1a in facies character and clast composition, being dominated by lava-rich ash-tuff lithofacies (Figs. 3, 4, 9a, and 12e). Its interpolated distribution shows deposits spread chiefly northward and southward, with maximum thicknesses in northern ravines (Fig. 9a). This pattern is consistent with dilute PDC deposition, and may record either BAF-related PDCs, with poorly preserved channel-confined facies (Macorps et al. 2018; Charbonnier et al. 2023; Capra et al. 2018, 2024), or another form of vent/conduit-clearing activity, such as hydromagmatic explosions (Suhendro et al. 2026). Layer 2 marks the transition to the first sustained fall phase. It remains rich in lava clasts and contains the highest lithic contents, but also includes significant scoria, indicating a stronger juvenile contribution (Figs. 2b, 3, 4, and 12e). Layer 2 preserves mainly massive, clast-supported facies and lacks evidence for major depositional interruptions. Its dispersal axis is directed northeastward (650 km2), unlike the more easterly dispersal of subsequent columns, suggesting that it records a distinct early column rather than the initial stage of layer 4 (Figs. 9d and 10a and b). Column-height estimates of 21 ± 3.8 km and MER values of 3.3 × 107 kg s−1 indicate a sustained sub-Plinian column (Figs. 11 and 12e).

The climactic stage began with deposition of layer 3, a matrix-supported tuff dominated by scoria and characterized by massive to diffuse and cross-bedded lithofacies, indicating emplacement from a dilute PDC (Figs. 3, 4, 9b, and 12f) (Branney and Kokelaar 2002; Kim et al. 2014). Although interpolated thicknesses are modest, layer 3 covers one of the broadest areas of any PDC in either eruption (221 km2), suggesting a highly mobile, unconfined current (Torres-Orozco et al. 2017a, 2018). Layer 3 may represent final cleansing preceding full establishment of a well-connected conduit, followed by deposition of layer 4, the largest and most widely dispersed fall deposit of the entire sequence (1500 km2), strongly dominated by scoria and defining a clear easterly dispersal axis (Figs. 2b, 3, 4, 9e, and 10c and d). An average column height of 29 ± 0.9 km, MER of 3.5 × 108 kg s−1, and magnitude of 4.6 identify layer 4 as a sustained Plinian column (Figs. 11 and 12f).

Following this climax, the eruption became progressively unstable, possibly due to advanced outgassing and waning gas-rich magma supply, as suggested by increasing dense lava clast contents (Figs. 2b and 12g) (Cioni et al. 2003, 2008; Shea et al. 2012; Torres-Orozco et al. 2025). Layer 5 records the most violent PDC phase and forms the most extensive PDC deposit (226 km2) (Fig. 9c). Despite its broad distribution, interpolated thicknesses reach only maximum 20 cm in restricted eastern sectors (0.03 km3 bulk volume). Its deposits are dominated by tuff lithofacies with the greatest internal variability of any layer, including massive, diffuse and cross-bedding, and locally clast-supported lenses, indicating a highly dynamic current (Figs. 3 and 4) (Belousov et al. 2007; Kim et al. 2014; Torres-Orozco et al. 2018). Layer 5 is followed by layer 6, another fall deposit dispersed eastward (430 km2) (Figs. 9f and 10e and f). An average column height of 24 ± 0.4 km and MER of 1.1 × 108 kg s−1 indicate that layer 6 represents a waning sub-Plinian column, most likely the declining remnant of the layer 4 Plinian phase after partial column collapse generated layer 5 (Figs. 11 and 12g). The eruption ended with deposition of layer 7, a concentrated PDC dominated by massive tuff, whose facies and geometry resemble the final collapse deposit of Lower layer 5 (Figs. 3, 4, and 12h).

Relationship between magmatic processes and eruptive behavior

The Lower Citlaltépetl deposits record interaction between compositionally distinct magmas. Juvenile products define a wide whole-rock compositional range, from dominating basaltic-andesite banded scoria and lava clasts to high-Si andesite pumice (Fig. 7a–l), showing clear mingling/mixing between these end members (Fig. 7) (Cas et al. 2023; Arce et al. 2026). Petrographically, scoria and pumice domains in banded clasts differ sharply: scoria is highly vesicular, dominated by spherical microvesicles, and containing amphibole + clinopyroxene + plagioclase, whereas pumice is overall vesicular, dominated by large coalesced vesicles, and containing plagioclase + orthopyroxene; both comprise abundant plagioclase sieve textures and zoning, indicating recharge and disequilibrium (Fig. 5h and i). Although abundant, amphibole fractionation does not appear to control the main bulk trends; instead, the evidence suggests that melt extraction, recharge, and magma mixing were more important. These relations suggest that the silicic component was a differentiated melt, likely extracted from a shallow evolved mush (Fig. 12b–d) (Sparks et al. 2022; Torres-Orozco et al. 2024), followed by increasing interaction with an ascending mafic magma that produced the dominant banded scoria during the climactic phase. Therefore, magma mixing was the main process driving the Lower Citlaltépetl sub-Plinian stage. The marked decline in banded scoria in layer 5, coupled with the greater abundance of individual pumice and scoria clasts, indicates that mixing diminished during the waning phase. This stage, characterized by reduced energy and column collapse, but lacking any increase in dense, outgassed lava fragments, likely reflects a shift from mingled to more homogeneous magmas, and progressive exhaustion of the evolved melt-rich domain (Cassidy et al. 2018; Sparks et al. 2022; Torres-Orozco et al. 2023).

By contrast, the Upper Citlaltépetl was fed by a compositionally homogeneous mafic-andesitic magma, with limited evidence of interaction with silicic melt. Bulk compositions cluster tightly within the andesitic field and show a narrow range (Fig. 7a–l). Trace-element compositions support this relative homogeneity: Upper Citlaltépetl forms a restricted cluster with limited dispersion in most ratios (Fig. 7d–k). Most Upper juvenile products are scoria; pumice is absent as a clear juvenile component, and most banded clasts are produced by internal textural heterogeneity within the scoria itself, ranging from microvesicular to microcrystalline to fibrous, vesicle-sheared scoria (Fig. 6c–f) (Shea et al. 2010, 2012; Torres-Orozco et al. 2023, 2024). Minor pumice-bearing banded clasts in the climactic phase may reflect faint silicic mingling, but petrography suggests that pumice could be rather accessory. Common clinopyroxene aggregates and rare olivine are consistent with a hotter, less evolved, and likely deeper mafic source (Fig. 7g). Thus, the main driver of the Upper climactic phase may have been strong rheological contrasts within the mafic magma from the onset of eruption (Fig. 12e–h), controlled by rapid magma ascent and decompression, strong bubble and melt coupling, relatively high crystallinity, and magma shearing (La Spina et al. 2021; Torres-Orozco et al. 2023; Bamber et al. 2024; Arzilli et al. 2026). Large-zoned and sieve-textured plagioclase record recharge and disequilibrium, whereas acicular to tabular microlites and microcrysts record rapid ascent and decompression (Arzilli et al. 2019). These conditions may explain the progression from early sub-Plinian to the climactic Plinian phase, and then to the increasingly unstable phases of layers 5–7 as gas supply waned and outgassing advanced.

Previous studies have proposed that effusive eruptions are favored by efficient magma mixing, which enhances heat transfer throughout much of a silicic magma body (e.g., Ruprecht and Bachmann 2010). This process raises pre-eruptive temperatures and reduces magma viscosity. In contrast, explosive eruptions may be associated with mafic recharge coupled with limited mingling and inefficient mixing, conditions that fail to supply sufficient heat to significantly lower the viscosity of the silicic magma (Ruprecht and Bachmann 2010). Our results from the Lower Citlaltépetl support this framework, indicating that highly explosive activity can be triggered by the injection of mafic magma combined with restricted mingling/mixing. Similar, relatively inefficient andesite-dacite magma mixing, producing marked banded scoria clasts, controlled the ~5 ka Jamapa pyroclastic flow at Pico de Orizaba (Arce et al. 2026). Furthermore, the transition from sub-Plinian to less explosive activity may reflect a shift from mingled magmas toward more compositionally homogeneous melts. This evolution implies a corresponding rheological transition that likely exerts a primary control on changes in eruptive style. These findings highlight a relatively underexplored aspect of magmatic systems: the role of rheological transitions between poorly mixed (mingled) and more homogeneous magmas, and how such transitions govern shifts in eruptive behavior.

Hazard implications

More than 500,000 people live within 30 km of Pico de Orizaba. The principal urban centers alone: Orizaba metropolitan area (~500,000 inhabitants), Coscomatepec (~60,000), Tlachichuca (~32,000), and Atzitzintla (~10,000), total approximately 602,000 people, excluding dispersed rural settlements (INEGI 2020). Several agricultural communities lie only 5–6 km from the summit (Fig. 1b), on older fall and PDC deposits, with limited services and poor evacuation roads. The Orizaba urban–industrial corridor, one of Mexico’s major manufacturing areas, lies 20–30 km southeast of the crater. Renewed activity comparable to the Lower or Upper Citlaltépetl eruptions could therefore affect both vulnerable proximal communities and densely populated distal sectors. Our results show that Pico de Orizaba can generate a broader range of explosive hazards than represented in previous studies and the existing hazard map (Sheridan et al. 2001). For the first time, we provide detailed distributions of several PDCs together with multiple sub-Plinian and Plinian fallout scenarios, showing that hazard varies with eruptive style, magnitude, column stability, current type, and plumbing-system processes. Even the smaller mapped PDCs (< 0.01 km3) are relevant because no previous study documented their distribution sufficiently for realistic scenario development.

Fall deposits record contrasting thicknesses, dispersal axes, areas (400–1500 km2), and reaches (Figs. 8 and 10). Upper layer 2 dispersed northeastward, Upper layers 4 and 6 eastward, and Lower layers 2+4 southeastward. Upper layer 4 represents the largest scenario, producing the thickest and most extensive deposit and a 30-km-high Plinian column; the remaining fallout units define smaller sub-Plinian scenarios (Figs. 8, 9, 10, and 11). In Fig. 11b–d, Lower layers 2+4 and Upper layers 2 and 6 overlap the broad range of sub-Plinian to lower-Plinian eruptions, whereas Upper layer 4 approaches large historical events such as the 1947 Hekla, 1913 Colima, and 1982 Chichón eruptions (Carey and Sigurdsson 1986; Mastin et al. 2009; Saucedo et al. 2010). Similar future activity could heavily affect the eastern and southeastern metropolitan and industrial sectors, while shifting winds could expose settlements in other directions. Fallout zonation should therefore distinguish smaller sub-Plinian from larger Plinian scenarios using different column heights, dispersal directions, areal coverages, and reaches. Comparable eruptive styles may also arise from contrasting plumbing conditions, from mafic-silicic mixing to decompression of rheologically heterogeneous mafic magma, and should not be represented by a single template.

The PDC deposits define equally contrasting hazards, including dome-collapse BAFs, dilute currents ranging from small, restricted to broad, highly mobile, and concentrated column-collapse currents (Figs. 8, 9, 10, and 12). Upper layers 3 and 5 demonstrate variable PDC energy and mobility, whereas Lower and Upper layer 5 show that PDCs of similar bulk volume (0.03–0.05 km3; Online Resource 1) may produce markedly different footprints. Lower layer 5 was confined and concentrated; Upper layer 5 was mainly unconfined, dilute, and more extensive. PDC hazard therefore depends not only on erupted volume and MER, but also on current behavior, mobility, topographic confinement, and generation mechanism. Both eruptions further show that PDCs may occur during vent clearing or dome collapse, waxing and unstable phases, partial column collapse, and final total collapse (e.g., Belousov et al. 2007; Shea et al. 2012; Charbonnier et al. 2023; Capra et al. 2024; Carazzo et al. 2025). Proximal settlements only 5–6 km from the crater and communities along valleys and ravines would be particularly difficult to evacuate if such currents accompanied renewed activity.

Future hazard assessment should therefore replace generalized envelopes with style-specific, multi-phase scenarios integrating fallout, contrasting PDC types, dispersal, transport mechanisms, topographic controls, and plumbing conditions. Such reconstructions are essential for land-use planning, evacuation design, and reducing loss of life around Pico de Orizaba.

Conclusions

This study shows that the Lower and Upper Citlaltépetl were two distinct early Holocene multi-phase andesitic explosive eruptions, composed of alternating fallout and PDC layers, which differed strongly in scale, dispersal, componentry, petrography, geochemistry, and eruptive behavior. The Lower Citlaltépetl was a smaller sub-Plinian eruption, with a total volume of 0.1 km3 DRE and magnitude 4.4, whereas the Upper Citlaltépetl was larger and more intense, with a total volume of 0.3 km3 DRE and magnitude 5. Individual fallout phases ranged from sub-Plinian to Plinian, with average column heights of ~21–30 km, and climactic Upper layer 4 representing the largest and most intense event recognized in this work.

The Lower Citlaltépetl is distinguished by abundant banded scoria and pumice, wider chemical compositional variability, and petrographic evidence for interaction between contrasting juvenile domains. These features suggest that its climactic phase was mainly driven by mingling/mixing between a shallow evolved high-Si andesite melt and an ascending mafic magma. In contrast, the Upper Citlaltépetl is dominated by scoria, has a narrower andesitic compositional range, and shows banded clasts produced mainly by internal scoria heterogeneity rather than magma mixing. Its eruptive behavior is therefore better explained by a larger, dominantly mafic-andesitic system in which gas-rich magma supply, rheological heterogeneity, strong bubble-melt coupling, rapid magma decompression and ascent, and magma shearing played the leading role. Together, both eruptions define two different pathways toward sub-Plinian and Plinian activity at Pico de Orizaba.

A major outcome of this work is the first detailed mapping of multiple PDC distributions at Pico. These maps show that PDCs ranged from BAF-related dilute PDCs to highly mobile broad currents, and concentrated column-collapse currents, each with different footprints. Likewise, fallout dispersal varied strongly between eruptions and even between layers, with northeastward, eastward, and southeastward axes and areal extents of ~400–1500 km2. These results show that Pico de Orizaba’s explosive hazard is more diverse than previously represented, and that future hazard assessment should move from generalized envelopes toward style-specific, multi-phase scenarios that include both fallout and contrasting PDC types.