Introduction

Persistent activity at open-vent volcanoes can provide excellent field laboratories to study the dynamics and evolution of different processes such as (i) the shallow conduit dynamics spanning from mafic lava lake-bearing to more evolved dome-forming systems (e.g. Harris 2008; Miwa et al. 2013; Scharff et al. 2014), (ii) the degassing and outgassing balance and its implications for understanding and modelling the magma feeding system (e.g. Aiuppa et al. 2010; Beckett et al. 2014; Edmonds et al. 2022) and (iii) transitions and controls on explosive–effusive volcanic eruption styles (e.g. Woods and Koyaguchi 1994; Ripepe et al. 2005; Cassidy et al. 2018). Open-vent volcanoes also allow calibration and association with characteristic geophysical signals allowing a better understanding of the conduit and source dynamics of explosive and effusive eruptions or modelling the plumbing systems (e.g. Ripepe et al. 2017; Magee et al. 2018; Gottschämmer et al. 2021). Such volcanoes are typically basaltic like Erta Ale, Etna or Stromboli (e.g. Allard et al. 1994; Clocchiatti et al. 2004; Sawyer et al. 2008), but can also be andesitic like Sakurajima or Colima to dacitic like Santiaguito (e.g. Mora et al. 2002; Harris et al. 2003; Miwa et al. 2013). Here, we take the period of persistent open-vent activity that occurred at Arenal between 1968 and 2010 to illustrate (i) the activity evolution and dynamics, (ii) the associated geophysical signals and (iii) the insights regarding the eruptive process behind the observed signals, while (iv) reviewing the literature for this long-lived and well-documented phase of open-vent activity at a basaltic-andesitic system.

On the morning of 29 July 1968, a new eruptive phase (VEI = 3) began at Arenal volcano, after several centuries of dormancy (Melson and Sáenz 1973; Borgia et al. 1988; Soto and Alvarado 2006). Although there are doubtful reports of possible minor eruptive activity in 1915 and 1922 (Soto et al. 1996), it is clear that the previous major explosive eruptive period occurred around 1430–1460 AD and that an effusive stage occurred around the same time. The volcano may even have remained active for several decades, or even into the seventeenth century (Soto and Alvarado 2006). Arenal has had highly explosive eruptions (plinian, subplinian and violent strombolian) every 200–600 years, alternating with minor eruptive phases, periods of inactivity and voluminous cone-building effusive phases (Soto and Alvarado 2006). The maximum Holocene eruption (VEI = 4) was subplinian to plinian type (Soto and Alvarado 2006).

Arenal is the volcano that has caused the most deaths in Costa Rica (82 in total, Alvarado et al. 2007) and the only one that has emitted lava flows in historical times. According to Soto and Sjöbohm (2005), the expected short-term hazards are tephra fall, pyroclastic density currents and small volcanic debris avalanches. In contrast, the long-term hazards are plinian eruptions, voluminous lava fields, mid- to large-volume volcanic debris avalanches and seiches in Arenal lake. A seiche is a standing wave in closed or partially closed body of water of finite length (Wilson 1972). This hazard was studied in Arenal by Hidalgo (1997) by modelling a pyroclastic density current or major rockslide entering the lake.

The open-vent activity at Arenal volcano (1968–2010) was particularly important for Costa Rica and thereby worthy of review, as the development of the capacity to measure, understand, track and monitor eruptive events led to the creation and evolution of a volcano observatory system. In 1974, the Instituto Costarricense de Electricidad (ICE) began monitoring Arenal and created an observatory that today is called the Observatorio Sismológico y Vulcanológico del Arenal y Miravalles (OSIVAM). During the 1980s, two other volcano-monitoring entities emerged. First, in 1983, the Red Sismológica Nacional (RSN) was founded following an agreement between Instituto Costarricense de Electricidad and the Universidad de Costa Rica (UCR). It then became a research programme within the UCR. In 1986, the Observatorio Sismológico y Vulcanológico de Costa Rica (OVSICORI) was founded at the Universidad Nacional Autónoma de Costa Rica (UNA). During the nearly 42 years of the eruption, Arenal became a natural laboratory for the study of volcanic processes in open-vent systems, making it the most studied volcano with the most complete eruptive record in Costa Rica.

The aim of this work is, therefore, to perform a retrospective analysis of the evolution of Arenal eruption (July 1968–October 2010) from a volcanological and geophysical point of view. We take advantage of the wide spectrum of geophysical, geological and geochemical studies carried out at Arenal volcano over the last 50 years, and in particular our own direct observations and field work carried out over the period 1990–2010. Information from local reports and publications was compiled and then integrated with unpublished geological field observations and geophysical data (seismological and acoustic) to allow a comprehensive view of the eruption. The new material presented and discussed here includes a summary of new geological observations on the emplacement of lava flows, the collapse of lava flow fronts and related block aprons and pyroclastic density currents, the provision of a revised map of the extent of the 1968–2010 deposits, a recalculation of the erupted volume and a reinterpretation of the evolution of the eruption. We also provide a preliminary analysis of volcano-tectonic seismicity and a recalculation of the seismic site effects. Finally, we provide an analysis of acoustic data obtained in 2005 that leads to new insights into the explosive mechanisms operating at Arenal. Our review is fundamental for proper hazard assessment and risk management at Arenal and may be useful for monitoring and understanding similar open-vent systems.

Geological setting

The Quaternary volcanic front of Costa Rica trends from north-west to south-east, 150-km inland of the Mid-American Trench. Subduction beneath Costa Rica is highly complex due to variations of the Cocos slab angle, morphology and petrological characteristics, as well as its interaction with the Caribbean plate and Panama microplate (e.g. von Huene et al. 1995; Barckhausen et al. 2001; Gazel et al. 2021). The Costa Rica volcanic front overlies a Mesozoic oceanic basement of Pacific and Caribbean plate origin and a Cenozoic volcano-sedimentary cover (Gazel et al. 2021). It is partitioned into three segments (Fig. 1a): the Guanacaste volcanic range at the north-western end, the Tilarán volcanic range and the Central volcanic range at the south-eastern end. The youngest volcanoes along the Tilarán volcanic range are Los Perdidos andesitic dome complex (86–83 ka), Chato volcano (48–3.7 ka) and Arenal volcano (7 ka to the present). They are aligned along a zone of crustal weakness-oriented NNW–SSE (Gillot et al. 1994; Soto and Alvarado 2006; Alvarado and Gans 2012).

Fig. 1
figure 1

a Simplified tectonic framework of Costa Rica and location of Arenal volcano. Abbreviations are: Guanacaste volcanic range (GVR), Tilarán Volcanic Range (TVR), Central Volcanic Range (CVR), the Central Costa Rica Deformed Belt (CCRDF) according to Denyer et al. (2013), Panama Fracture Zone (PFZ) and Nazca plate (Np). Black continuous and dotted lines represent faults. Triangles denote reverse faulting and points towards the overriding plate. Bathymetry and topography from Ryan et al. (2009). b Arenal volcano topographic surface and main features: in circles the old crater (D) and craters generated during the onset of the eruption, 29–31 July 1968 (A, B, C) and extension of the resulting 1968–2010 lava field (blue contour). Main labelled faults (red colour) are: Aguacaliente River fault (ARF); Danta fault (DF). The dotted lines represent the fault traces covered by the lava field. Instituto Geográfico Nacional provided the Arenal digital elevation model based on 2014–2019 images. Local coordinates EPSG 5456 Lambert North

Full size image

Arenal is located at the NW end of the Tilarán volcanic range, within a ~ 70-km zone of decreased Holocene volcanic activity that coincides with a change in the morphology and subducting angle of the Cocos plate (Fig. 1a). The subducting slab has a smooth surface and steeper angle off the north Pacific coast and a seamount domain off the central Pacific coast with an oceanic crust overprinted by Galápagos tracks (e.g. Stoiber and Carr 1973; Wade et al. 2006; Dinc et al. 2010; Gazel et al. 2021). Different studies agree that crustal thickness beneath Arenal ranges from 35 to 40 km (Matumoto et al. 1977; MacKenzie et al. 2008; Lücke 2014).

Arenal was formed within a tectonic depression with a Mio-Pleistocene volcanic basement, filled with lacustrine and alluvial deposits and regional tuffs, with an age of ~ 20 ka possibly from the Chato volcano (Soto and Alvarado 2006). Arenal's edifice construction started with the deposits of a phreatomagmatic eruption (AR-1) around 7 ka, which were emplaced on a well-developed lateritic paleosol (Soto and Alvarado 2006). This was followed by 21 explosive eruptions of VEI ≥ 3 (AR-2 to AR-22), and the emplacement of at least five lava flow fields up to 1968, building the present cone (Borgia et al. 1988; Soto and Alvarado 2006).

Review of the eruption and associated activity

The July 1968- October 2010 eruption: timeline

Arenal's 42-year-long eruptive phase was part of six long-lasting (> 30 years) eruptions in Central America spanning the twentieth and twenty-first centuries, together with Santiaguito, Fuego, Pacaya, Izalco and San Miguel volcanoes (Alvarado et al. 2007). This long-lasting activity began on 29 July 1968, when Arenal awoke violently from its dormancy with an explosive eruption that lasted three days in an event characterised by an initial lateral blast followed by pyroclastic density currents, ashfall and ballistics. The opening phases of the eruption have been widely documented by Melson and Sáenz (1968, 1973) and Minakami et al. (1969), as well as by Alvarado et al. (2006), with the eruption itself continuing until October 2010 and following seven stages. The location of the vents along a radial fissure is consistent with the west-directed blast during the onset of the eruption (Melson and Sáenz 1968; Melson and Sáenz 1973; Alvarado et al. 2006) and the direction of maximum horizontal compressive stress (López 1999; Alvarado et al. 2006).

Stage 1 July–September 1968

A lateral blast opened three new vents, roughly aligned east–west (Figs. 1b and 2a), and called A (~ 1060 m a.s.l.), B (~ 1160 m a.s.l.) and C (~ 1400 m a.s.l.) (Melson and Sáenz 1968, 1973; Minakami et al. 1969). It is not possible to determine whether their formation was simultaneous or successive to the initial blast (Alvarado et al. 2006). The block ejection reached over 5 km away from the lower crater A (Melson and Sáenz 1968). The pre-existing crater was called crater D (Figs. 1b and 2a). The violent explosive phase was followed by small explosions and declining fumarolic activity through August and mid-September, until stronger explosions occurred between 13 and 18 September (Alvarado et al. 2006).

Fig. 2
figure 2

adapted from Fig. 2 of Almendros et al. 2014) and d 2000–2010 (adapted from Fig. 8b of Valade et al. 2012). Dark grey represents the pre-existing cone and light grey is the new cone formed since 1974. Yellow to orange colours represent the volcanic conduits and brown colour colder lavas

Schematic E–W cross sections showing four main stages of Arenal evolution: a 1968, b 1974, c 1987–2000 (

Full size image

Stage 2 September 1968-August 1973

The opening phase was followed by a five-year-long period of effusion of basaltic–andesitic lavas from crater A. This proceeded at a typical dense rock equivalent discharge rate of 0.93 m3 s−1 (Wadge et al. 2006). In this review, we define discharge rate as the volume flux averaged over a given time period (Harris et al. 2007).

Stage 3 August 1973-March 1974

This period marked a seven-month hiatus during which almost no effusion of lava occurred. Fumarolic activity did, though, continue, especially in crater C from which lava flow effusion resumed in March 1974. This marked the onset of Stage 4 (Fig. 2b).

Stage 4 March 1974-mid-1984

As part of Stage 4, the focus of eruptive activity shifted from crater A to crater C (Wadge 1983). Blocky lava flows were almost continuously emitted. An explosive event generated a pyroclastic density current in June 1975 (Matumoto and Umaña 1976; Van der Bilt et al. 1976; Alvarado and Soto 2002) and small explosive phases occurred in August 1980 (Güendel and Malavassi 1980) and May–June 1981 (Cheminée et al. 1981). The dense rock equivalent discharge rate for the first three years was as for Stage 2, i.e. 0.93 m3 s−1, but they dropped in 1977 to 0.25 m3 s−1 and remained at these levels for the remainder of Stage 4 (Wadge et al. 2006).

Stage 5 mid-1984–1987

In mid-1984, Arenal entered an explosive phase. This was accompanied by an increase in the dense rock equivalent discharge rate to 0.61 m3 s−1 (Wadge et al. 2006). Continuous effusion of lava flows continued into 1987 (Fig. 3a and b).

Fig. 3
figure 3

a Typical lava flows from cone C on the W side, with well-developed channels, lateral levees and block-rich fronts, dividing in two or more lobes. Photograph from 29 October 2007. b Arenal from the W. In the foreground, the deposits of the block-and-ash flow of 30 July 1968. In the middle, blocky lava flows from the 1980s and 1990s, with different vegetation covers, according to age. In the background, the cone C, with blocky lava flows (dark grey) and block aprons (light grey). Photograph from 22 June 2009. c Blocky lava flow from 1995 on the W side of Arenal in the foreground, 1968 pyroclastic density currents deposits and recent epiclastics with gullies in the middle, and the Arenal reservoir in the background. Photograph from 6 March 2007. d Lava flows, block aprons and gullies developed on the N side of Arenal, covering the dense forest. Photograph from 27 February 2008. e Arenal seen from the NNW. On the far left, part of Chato volcano. Cone D is on the left and cone C on the right. The lava flows and block aprons are in grey colour. Photograph at the very end of the eruption, on 7 October 2010. f Arenal degassing seen from the S. Cone D on the right, covered by dense vegetation. Photograph from 30 January 2008. g Pyroclastic density currents developed from a lava front collapse towards the SW flank of Arenal on 1 December 2008. Photograph courtesy of Luis A. Madrigal (Instituto Costarricense de Electricidad). h Close-up of the top of the SW side of cone C. Vigorous degassing from the lava pond, very short and slow lava flows on a steep slope (so-called dome-like), and development of block aprons by collapse of the lava fronts. Photograph from 6 March 2007. i One of the very last ash and gas explosions with weak and short ash columns. Photograph taken at 16:39 on 7 October 2010, from the E, near La Fortuna town. All photographs but g were taken by Gerardo J. Soto

Full size image

Stage 6 1987–2000

Stage 6 was marked by the accumulation of lava inside crater C which, although becoming capped by a thick stable crust, fed longer reaching lava flows (≥ 0.7 km-long) (Fig. 2b and Fig. 3b and c). The longest lava flow reached 3.2 km in 1992. Stage 6 was also marked by major pyroclastic density currents (2–3.2 km long) associated with crater wall collapse (Alvarado and Soto 2002), intense explosive activity and the formation of strombolian “hornitos” associated with short-lived low lava fountains in August 1996 (Soto and Arias 1998), April 1999 and September 1999 (Arroyo et al. 1999). The dense rock equivalent discharge rate ranged between 0.24 and 0.31 m3 s−1 (Soto and Arias 1998; Wadge et al. 2006). During the 1980s and 1990s, the lava flows were continuous bodies of lava enclosed by debris or levées and frontal zones composed of a distal accumulation of blocks and a protrusion of lava that emerged from the debris (Linneman and Borgia 1993). The vertical section of the lavas consisted in a bottom of debris followed by an unvesiculated core and a top of debris (Cigolini et al. 1984).

Stage 7 2000-October 2010

Stage 7 was marked by the effusion of shorter (< 0.6 km long) lava flows and the building of a complex structure over crater C (Fig. 2d). After 2000, the lava discharge rate dropped to 0.1 m3 s−1 and remained at such low levels until October 2010. Stage 7 was also associated with reduced explosive activity, formation of small dome-like structures and occasional strombolian “hornitos” (Fig. 3d-i).

During the last stage of the eruption, a decrease in seismic activity was also observed. The real-time seismic energy measurement (RSEM) gradually decreased until July 2010 when it dropped to background noise levels (Fig. 4a). Only a few long-period events were observed after October 2010, and the last harmonic tremor was recorded on 30 October 2010 at 15:14 UTC. Furthermore, the dominant frequency remained stable around 3 Hz until March 2010 when it gradually decreased up to 2 Hz (Fig. 4b). The last explosions were witnessed on 7 October 2010 (Fig. 3i) and lava effusion stopped by late October. Thereafter, the activity consisted of falling blocks from the inactive lava flow front until early 2011 (Fig. 3e).

Fig. 4
figure 4

Seismic activity recorded at ARE1 station from Red Sismológica Nacional (located 1.6 km westward from the active crater) from June 2005 to December 2010. a Real-time seismic energy measurement (RSEM, De la Cruz Reina et al. 2001). b Dominant and maximum frequencies following Carniel and Iacop (1996). All estimates were carried out within the range 0.5–10.5 Hz and using a 60-s sliding window. Results were smoothed using a moving average to obtain a point every 15 days. Shaded grey vertical lines enclose periods with no data due to instrument failures. In June 2006, the seismometer was changed from a Lennartz LE-3Dlite to a Guralp 6TD 30 s

Full size image

Evolution of crater C

Morphology

During the onset of the eruption in July 1968 (Stage 1), three craters (Figs. 1b and 2a) were formed. Crater A remained active (with lava effusion until August 1973 and then only fumarolic) until March 1974, when activity migrated to crater C, which was about 400 m higher than crater A (Stage 4). Crater A was later covered by the lava flows emitted from crater C. While the morphology of crater A is not known in detail, crater C and the cone it grew on were described more extensively throughout the eruption. Melson and Sáenz (1968) and Melson and Sáenz (1973) roughly drew both craters having an elliptical geometry oriented N 80° W. The major and minor axes of crater A were ~ 400 and ~ 240 m, respectively, whereas the major and minor axes of crater C were ~ 160 and ~ 80 m, respectively (Melson and Sáenz 1968; Melson and Sáenz 1973).

The almost continuous emission of lava flows from crater C since 1974 progressively built a parasitic cone C, which grew from ~ 1400 m a.s.l. up to 1580 m a.s.l. in 1986 (Fig. 5a). There are no much descriptions of this growth, but for crater shape and dimensions in the early 1980s. The cone grew faster once the explosive activity resumed in 1984, making it wider and higher, and transforming it from a parasitic cone to a twin-cone composed by lava flows and pyroclastics, and then, transforming Arenal from a simple-cone edifice to a sub-cone edifice (Fig. 3e-f) (cf. Grosse et al. 2009).

Fig. 5
figure 5

Eruptive activity of Arenal in 1986–2010. a Discharge rate (black) based on data obtained from Soto and Arias (1998), Wadge et al. (2006) and estimations by one of the authors (GJS) for the period 2005–2010. Crater C elevation changes (red) from Wadge et al. (2006) and normalised (blue). b Number of explosions per hour (bars) recorded at the reference station FOR (located 3.5 km E from the active crater) from Observatorio Sismológico de Arenal y Miravalles and ash fall collected at 1.8 km from the active crater (red line) from 1992 to 2006 (data from Soto 1998) and Vargas et al. 2006). c Number of pyroclastic density currents per year (bars) and maximum length (blue circles and lines). d Documented descending flanks of the pyroclastic density currents

Full size image

In 1982, crater C had a 60-m-diameter circular shape and was filled with lava (Cigolini et al. 1984). By 1988, it doubled in diameter (~ 120 m), while the extension of the summit complex was ~ 250 m north–south and ~ 170 m from east to west (Cole et al. 2005). Alvarado and Soto (2002) suggested that crater C was funnel-shaped to explain the crater widening, but the dimension of the vent at the floor of the crater is not known. The collapse of the crater rim that occurred in 1993 led to a major pyroclastic density current on 28 August and left a horseshoe-shaped opening about 80 m deep (Alvarado and Soto 2002). In 1998, the C cone was ~ 250 m north–south and ~ 170 m east–west, and the crater was ~ 120 m in diameter (Cole et al. 2005). In March 2000, crater C was ~ 150 m diameter and had a shallow amphitheatre morphology open to the west (Cole et al. 2005). Between 2000 and 2010, the development of the active cone was complex, being dominated by an accelerated growth (mainly during 2002–2005, Fig. 5a). In May 2016, the summit crater measured ~ 100 m north–south and 150 m east–west (based on Google Earth imagery).

Wadge et al. (2006) indicated that crater C reached a final elevation of 1715 m in 2005 (Fig. 5a). However, the topographic map provided by the Costa Rican agency for mapping Instituto Geográfico Nacional (IGN), based on images taken between 2014 and 2019, reveals that cone C reached a final elevation of 1670 m a.s.l. This means that the total growth of the new cone was 270 m since 1974. Furthermore, considering that crater C and crater D were almost at the same level in 1989 (Linneman and Borgia 1993), it is possible to normalise the data of Wadge et al. (2006). For this purpose, linear trends were obtained from the original data of Wadge et al. (2006) and from that between the elevation of cone C reported in 1989 by Linneman and Borgia (1993) and in 2014 by the National Geographic Institute. The difference between the two lines was then obtained and subtracted from the Wadge et al. (2006) data. We can observe that the growing trend of cone C is less accelerated than reported by Wadge et al. (2006) (Fig. 5a).

Dynamics

During the 1980s, crater C was described as a nearly circular pit filled with basaltic–andesitic lava (Cheminée et al. 1981; Global Volcanism Program 1981; Cigolini et al. 1984; Alvarado and Soto, 2002) from which pulsating glows were observed (Barquero, 1980). The lava overflowed from the crater rim (Fig. 2b) and moved down along the northern and south-western flanks of the volcano (Cigolini et al. 1984; Alvarado and Soto, 2002). Between June 1979 and April 1980, Cheminée et al. (1981) and Global Volcanism Program (1981) described that crater C was plugged with lava forming a dome-like structure from which lava effused and gases escaped through radial fissures. In March 1981, violent degassing of the upper part of the vent expelled large plumes of vapour, sometimes including tephra and bombs, which fell 100 m from the crater rim (Cheminée et al. 1981; Global Volcanism Program 1981). This activity was accompanied by loud detonations. Moreover, Cheminée et al. (1981) and Global Volcanism Program (1981) reported that, once the dome was deflated, the fissures sealed, building up pressure and generating small explosions ~ 0.2 km high. A similar behaviour was observed at Lascar volcano between 1984 and 1995 (e.g. Openheimer et al. 1993; Matthews et al. 1997; Wooster and Rothery 1997). Direct observations of the crater surface (Cheminée et al. 1981) described a uni-directional flow from the side of the crater with intense degassing, upwelling and bursting towards the opposite direction from where lava flowed out of the crater and downslope. This is similar to what has been observed at Halema'uma'u, Kīlauea (Patrick et al. 2019) and Villarrica (Lev et al. 2019). In addition, there were areas of profuse gas emission and bubble bursting (Cigolini et al. 1984). In 1983, another expedition to the crater C described gas emission along its edge (Benoit and McNutt 1997).

Between 1984 and 2000 (stages 5 and 6), the molten lava hosted inside crater C continued to form a cap which in turn generated frequent explosive activity (Melson, 1989; Cole et al., 2005; Alvarado 2011). During that period, the walls of crater C collapsed outward several times generating pyroclastic density currents and subsequent lava overflow (Alvarado and Soto 2002; Cole et al., 2005). After one of these episodes, in August 1993, the lava effusion point shifted from south to north (Wadge et al. 2006). In addition, multiple vents along the crater rim were commonly observed with different eruptive activities. In October 1994, there were two pyroclastic cones in the crater, one in the north from which lava effusion occurred, with profuse degassing, and another in the south that produced minor and moderate explosions and soundless ash emissions (Soto, 1997). In 1996, at least two vents were observed, one of which was the source of lava flow and over which a 30-m-high spatter cone was formed, and the other was the source of explosive activity generating plumes ~ 2 km high and ejecting blocks reaching 1 km high (Soto and Arias 1998). Occasionally, intense pulsating glows (every ~ 40 s) could be also observed over the crater without explosive activity (Soto and Arias 1998). According to our observations, between 2000 and 2010 (Stage 7) crater C became a growing summit complex characterised by multiple small vents, formation of dome-like structures and spatter cones (“hornitos”), punctuated by a mild explosive activity (Fig. 2d).

The hosted lava within crater C: a lava pond?

The molten basaltic–andesitic lava hosted inside crater C and located directly over the underlying conduit was commonly described as a “lava pool” (e.g. Alvarado and Soto 2002; Alvarado et al. 2003, 2010; Cole et al. 2005; Lesage et al. 2006; Bertolino et al. 2009). The term of lava pool or pond has been used in the context of basaltic volcanism to describe an accumulation of lava within a topographic depression (e.g. Tazieff 1994; Head and Wilson 1989; Vergniolle and Mangan 2000; Wolff and Sumner 2000). Lava ponds can be generated by lava fountains where lava is trapped at the vent by the walls of a developing spatter or cinder cone (e.g. Head and Wilson 1989; Vergniolle and Mangan 2000; Wolff and Sumner 2000), or located adjacent to, or down flow from, the source vent resulting from the lava flowing into a previously formed depression (e.g. Wilson and Parfitt 1993, Stovall et al. 2009; Harris 2008; Patrick et al. 2011). Such features are better known as inactive or passive lava lakes in which there is no convection and lava remains within the confining pit to stagnate and cool in place (e.g. Harris et al. 2005; Harris 2008; Lev et al. 2019). An example of an inactive lava lake was the ‘Alae pit crater in Kīlauea (Swanson et al. 1972). Therefore, it seems that lava pool or pond is not a proper term to describe the hosted lava within crater C. Several interpretations of this feature can be considered: 1) crater C could have hosted an active lava lake such as Masaya, Ambrym, or Erta Ale (e.g. Duffell et al. 2003; Harris et al 2005; Allard et al. 2016), or 2) it was just the uppermost part of a conduit as has been suggested for Erebus (Vergniolle and Bouche 2016; Vergniolle and Métrich, in revision, this issue) and Villarica (Witter et al. 2004), or 3) it was a flank feature sitting over an active vent and feeding further discharge such as the effusive vents observed during the 2002–2003 flank eruption at Stromboli volcano (Calvari et al. 2005). This is still an open question.

A convective regime allows a lava lake to remain active (e.g. Tazieff 1994; Harris 2008; Lev et al. 2019) and has been invoked to explain the unbalance between the large volumes of magma degassed and the small volume of magma erupted at persistently degassing volcanoes such as Izu-Oshima, Villarrica or Erta Ale (e.g. Kazahaya et al. 1994; Harris et al. 2005; Palma et al. 2011). Comparisons of SO2 petrological budget (based on melt inclusions) and SO2 flux (correlation spectrometry, COSPEC) were made at Arenal by William-Jones et al. (2001) and Wade et al. (2006). Williams-Jones et al. (2001) obtained a petrological emission of ~ 0.41 Mt based on clinopyroxene- and plagioclase-hosted melt inclusions which was 3–9 times lower than the 1.3 Mt in average obtained from COSPEC (and as high as 3.9 Mt, considering windspeed and other caveats). However, Wade et al. (2006) considered that the melt inclusions of Williams-Jones et al. (2001) were dacitic and hence considerably more degassed than the ones expected from a basaltic-andesitic eruption. Instead, Wade et al. (2006) considered their mafic to intermediate inclusions to be more representative of the modern Arenal magma and obtained SO2 estimates between 500 and 1500 ppm, which is within the range of those obtained with COSPEC by Williams-Jones et al. (2001). In general, this suggests that Arenal is being continuously supplied by fresh hot undegassed magma (Williams-Jones et al. 2001). Thus, there is probably no apparent unbalance between degassing and magma erupted. However, these results should be taken with caution as there is not enough data to characterise in more detail the SO2 flux variations throughout the entire eruptive stage of Arenal and further analysis should be carried out in this respect. Other observations like the uni-directional flow (Cigolini et al. 1984), glowing, “hornitos” and spatter cones are not conclusive of a convective regime as they are also observed during lava effusion (e.g. Head and Wilson 1989; Vergniolle and Mangan 2000; Calvari et al. 2005). While magma convection may operate in low-viscosity basaltic systems, this mechanism may not be important in more evolved volcanic systems, due to the extensive degassing-induced crystallisation in the conduit (Edmonds et al. 2022). Further research and modelling is needed to better constrain the properties and dynamics of the so called “lava pool” and the conduit at Arenal.

Explosive activity

The explosive activity during the 1968–2010 eruption is poorly understood in terms of mechanisms, products and temporal variations. The products of the 29–31 July 1968 explosive onset are described by Alvarado et al. (2006), while available information of the open-vent phase regarding eruptive column heights and ballistics can be found in Table S1 (Supplementary Material 1), and the few textural descriptions of explosion products can be found in Barquero and Alvarado (1989a), Dellino and Alvarado (1992) and Cole et al. (2005).

The main explosive phase that occurred between 29 and 31 July 1968, which is the onset of the 1968–2010 eruption, had a VEI of 3 and fed columns up to 10 km in height. This opening phase continued until mid-September, but the VEI dropped to 1–2 and lava effusion began (Alvarado et al. 2006). Thereafter, the eruption was mostly effusive until 1984, except for brief periods of explosive activity in 1975, 1980 and 1981. Frequent explosive activity started in mid-1984 and varied between one and six explosions per hour according to monitoring data collected between 1986 and 2004 (Fig. 5b). Cole et al. (2005) reported 250 gm−2 day−1 of tephra emitted in 1989, reflecting intense explosive activity that generated ash plumes reaching 1.5–2 km between 1984 and 1989 (Table S1, Supplementary Material 1). Tephra fall was collected by Observatorio Sismológico y Vulcanológico de Arenal y Miravalles at a site located 1.8 km west (735 m a.s.l.) of the active crater between 1992 and 2006 (Soto 1998; Vargas et al. 2006). This collection system showed a peak in August 1996 at 64 gm−2 day−1 and another in June–August 1997 at 54 gm−2 day−1. These peaks were associated with frequent and stronger explosions that produced 1–1.5-km-high ash-laden columns and ballistic blocks. The mass of tephra was < 10 gm−2 day−1 from 1998 to May 2006, when the last samples were obtained (Fig. 5b). From 2006 to 2010, the quantities recorded at the collecting site were only traces amounts. The ballistics ranged from 0.5 to 4 m in diameter and reached distances of 0.5 to 2.1 km from crater C between 1984 and 1992. These generated impact craters, with widths and depths of 0.8–2 m and 0.2–4 m, respectively (Table S1, Supplementary Material 1).

Melson (1989) and Barboza and Melson (1990) provided a detailed description of the explosive activity carried out from 1 to 13 April 1989, based on the types and intensity of the sounds and the associated products. On this basis, they distinguished three types of events: 1) explosions, 2) “whooshes” and 3) “chugs”. The explosions were the most intense events and had a short energy release (< 2 s). They were accompanied by blocks and bombs, 0.3–1-km-high plumes with low volumes of tephra and infrequent pyroclastic density currents. The strongest explosions occurred after a long quiescence time (~ 2 h). “Whooshes” were small sustained events that generated jet engine sounds with emission of blocks and bombs, also described as bomb fountain by Melson (1989), and fine tephra. “Chugs” were rhythmic gas emissions that sounded like a steam locomotive and with or without ejection of small amounts of tephra (Melson 1989; Barboza and Melson 1990; Benoit and McNutt 1997). “Chugging” has also been observed at other volcanoes such as Karymsky and Sangay (Johnson and Lees 2000), Reventador (Lees et al. 2008) and Fuego (Lyons et al. 2010). “Whooshes” frequently preceded the chugging, as has also been observed in the Karymsky, Sangay and Reventador (Johnson and Lees 2000; Lees et al. 2008).

The explosive activity at Arenal has been described as strombolian by most authors (Barquero et al. 1984; Barquero and Alvarado 1989a; Dellino and Alvarado 1992; Williams-Jones et al. 2001; Szramek et al. 2006). In contrast, few authors interpreted the explosive activity as vulcanian (Cole et al. 2005; Mora et al. 2013). Johnson (2004) and Alvarado (2011) considered that Arenal explosive activity could be both strombolian and vulcanian. To date, there is no consensus on whether the Arenal explosive activity was purely strombolian or vulcanian. We will address this aspect later in the discussion.

Pyroclastic density currents were also common at Arenal (Fig. 3g). For this review, we completed and revised the catalogue of pyroclastic density currents obtained by Alvarado and Arroyo (2000) by carefully reviewing the eruptive activity reports from the local observatory network and reports published by the Global Volcanism Program, as well as the studies of Alvarado and Soto (2002) and Cole et al. (2005). This updated catalogue is given as Table S2 in Supplementary Material 2. From 1986 to 1991, pyroclastic density currents occurred mainly by column collapse at a rate of < 20 events per year and with run-outs of 0.5–1.5 km towards the NW and S flanks (Fig. 5c and d). Between 1991 and 2001, the largest pyroclastic density currents (on 28 August 1993, 5 May 1998, 23 August 2000 and 24–26 August 2001) were generated by the collapse the crater C wall combined with the overflow of the lava pond (Alvarado and Soto 2002; Alvarado et al. 2003, 2010). The run-outs and velocities of these large pyroclastic density currents were 2–3.2 km and 11–33 m s−1, respectively, while subsequent lava flows reached velocities between 30 and 125 m day−1 and distances from 700 to 500 m the first day (Alvarado and Soto 2002).

Lava front collapses produced block-and-ash-type pyroclastic density currents during the last stage of the eruption (Stage 7: 2000-October 2010). Once the lava flow front collapsed, a depressurisation flow interior occurred, to generate the pyroclastic density current (e.g. Rose et al. 1976; Fink and Kiefer 1993; Harris et al. 2002). Consequently, the event changed from a rocky slide to a pyroclastic density current, over a very short period of time (less than one or two seconds). Similar mechanisms have been observed from the dacitic lava flow of Santiaguito in Guatemala (Rose et al. 1976; Harris et al. 2002; Rhodes et al. 2018) and the dacitic dome of Unzen in Japan (e.g. Suto et al. 1993; Yamamoto et al. 1993; Cordonnier et al. 2009). The pyroclastic density currents resulting from lava front collapses had run-outs < 1 km and were very frequent in the period 2001 to July 2010. From 2001 to 2005, pyroclastic density currents associated with flow front collapse flowed towards the N-NNE and NW, whereas from 2005 to 2010, the pyroclastic density currents flowed to the SW, S and SSE, depending on the direction of the lava flows.

Lava flow descriptions

The almost continuous effusion of lava from 1968 to 2010 at Arenal, resulted in a 4 km long and 4.5 km wide lava field that became part of the group of historical basaltic and basaltic–andesitic lava fields, along with other cases such as Kilauea and Mauna Loa, Etna and Piton de la Fournaise, among others (Kilburn 2000). This made Arenal a field laboratory for studying the morphology and dynamics of andesitic lava flows as it is evidenced by several studies such as Bennett and Raccichini (1977), Borgia and Cigolini (1980), Borgia et al. (1983), Cigolini et al. (1984), and Linneman and Borgia (1993). These studies provided fundamental information for modelling the morphology and dynamics of lava flows and lava flow fronts in terms of rheological and topographic constrains as well as understanding the developing of lava fields, particularly those of ‘a’a and blocky flows (e.g. Baloga and Pieri 1986; Dragoni et al. 1986, 2005; Kilburn and Lopes 1991).

Between 1990 and 2010, extensive mapping of the lava flows was done with binoculars and benchmarks previously set as references above the elevation of 900 m, which is the steepest and most dangerous area of the cone. Below this elevation, the mapping was done by basic topographic surveys. This allowed a complete description and characterisation of the lava flows, hence a better calculation of the discharge rates. The long lava flows from early 1990s (up to 3.2 km long) shifted to much shorter flows during the last stage of the eruption (2000-October 2010), when lavas were emitted at a rate of 0.1 m3 s−1 (a third of the rate in most of the 1990s) (Fig. 5a). The lava flows did not extend beyond the upper parts of the active crater (< 0.6 km) and moved along very steep slopes (> 35°), generating instabilities of the lava fronts that led to frequent incandescent rock slides, and eventual pyroclastic density currents, when the detached volume was high enough (usually > 104 m3). Such incandescent slides and the generation of pyroclastic density currents have been described elsewhere, as Stromboli (Lodato et al. 2007) and Santiaguito (Rose et al. 1976), respectively. The reworked and abraded rolling blocks (decimetric to metric) formed fans supported by grains or by sand. These deposits are very thick (tens of metres) in the middle part of the volcano, typically between 1 and 2 km from the active crater. The texture of these deposits is usually heterogeneous, facially and granulometrically, depending on the volumes detached from the flows.

Lava flow field and final volumes erupted

Wadge (1983) and Wadge et al. (2006) conducted estimates of the volume of the lava flow field and discharge rates at Arenal from topographic data collected by radar interferometry, stereo-photogrammetry and field survey. Wadge (1983) made the estimates for the period from 1968 to 1980, which were then recalculated and extended to 2004 by Wadge et al. (2006). Macfarlane et al. (2006) conducted short-term all-weather volcano topography imaging sensor (AVTIS) radar measurements in April and May 2005 to determine topographic changes associated with the advance of the lava flows and to estimate the discharge rate.

Although highly accurate, the topographic maps of the lava flow field from Wadge (1983) and Wadge et al. (2006) were obtained by different techniques and periods which makes estimates not necessarily uniform. For this review, we thus conducted a new estimation of the volume of the lava flow field and of the 1968–2010 erupted material by differencing the 1961 (1:50 000 scale map), 2005 (1:25 000) and 2014–2019 (1:5000) topographies provided by Instituto Geográfico Nacional. For all maps, we calculated 10-m-resolution elevation models by using triangulated irregular network (TIN) interpolation and then corrected by filling the no data areas by interpolating from valid pixels around the edges of the area. The three elevation models were referenced to the same CRTM05 coordinate system, which is the official one for Costa Rica. As no well-distributed reference points were available, the coordinates and elevations of the Instituto Costarricense de Electricidad dry tilt stations used during the 1980s and 1990s (Alvarado et al. 1988a; Soto 1991), as well as those of the seismic antennas deployed in 1997 (Métaxian et al. 2002; Mora 2003, Mora et al. 2004, 2006; Almendros et al. 2012, 2014), and 2008 (Guerrier 2012) were used to assess, approximately, the accuracy of the elevation models. The location and elevation of the dry tilt sites and seismic antennas were measured by the Topography Department of the Instituto Costarricense de Electricidad, so they are sufficiently precise for this purpose. A total of 131 points were obtained, of which 91 are on the western and south-western flanks and the rest are distributed on the northern and eastern flanks of the volcano (Table S3, Supplementary Material 3). We obtained a root mean square error (RMSE) (Pakoksung and Takagi 2016) values of 14.35, 7.3 and 8.15 m for the elevation models of 1961, 2005 and 2014, respectively, considering all the points (Table S4, Supplementary Material 3). If we use only the points situated on the northern (outside of the lava flow field) and eastern flanks, which are not strongly affected by the deposits of the eruptive activity, we obtain root mean square error values of 21.70, 3.10 and 3.60 m for the 1961, 2005 and 2014 elevation models, respectively (Table S4, Supplementary Material 3). In this last case, for the 2005 and 2014 elevation models the root mean square errors are improved, whereas it increases for the 1961 elevation model.

To estimate the volume of the lava flow field and the total volume of the 1968–2010 eruption products, we used the contours shown in figure S1 in Supplementary Material 3. To account for the boundary effect, we used a buffer of 20 m which represent uncertainties less than 1%. The results of the volume calculations are shown in Table 1. The correction of the volumes to dense rock equivalent was conducted considering global porosity of 19.6% by prorating 65% lavas at 14% porosity and 35% of pyroclasts, epiclasts and blocky fans at 30% porosity. We consider this proportion reasonable as the contribution of pyroclastic density currents and tephra fall is not quite minimal, but only into the lava field (Wadge et al. 2006). The porosity values for lavas and pyroclasts have been chosen following Wadge et al. (2006) and Alvarado et al. (2006). The dense rock equivalent volume of 527 ± 58 Mm3 for the lava flow field is consistent with the results obtained by Wadge et al. (2006). The additional volume added from 2005 to 2010 is about 25 ± 3 Mm3 which is in the order of the expected volume considering ~ 0.1 m3 s−1 (Wadge et al. 2006; Taylor and Soto 2010) for that period. The total dense rock equivalent volume change along the volcanic edifice is 732 ± 75 Mm3.

Table 1 Estimates of the lava flow field and total volume erupted at Arenal volcano
Full size table

If we compare the volcano surface before and after the eruption, we notice that most of the new volume is distributed towards the western and north-western flanks, where lava flows, pyroclastic density currents and ash fall accumulated (Fig. 6). To a lesser extent, there are subtle changes to the east and south-east flanks due to erosion and reworking of erupted material. The maximum elevation change is 292 m related to the maximum thickness of the lava flow field area corresponding to the cone C. In addition, the summit, which was almost perfectly conical before 1968, became asymmetrical due the construction of cone C and began to exceed the height of crater D in 1991 (Figs. 3e-f and 5d). While the elevation of crater C increased as described before, the elevation of crater D decreased by ~ 30 m, mostly due to erosion generated by explosive activity and bombardment of large ballistics from crater C, as it was witnessed by the authors through the eruption (i.e. Global Volcanism Program 1988) (Fig. 6).

Fig. 6
figure 6

Comparison of Arenal elevation models before the eruption and after the eruption. The following features are shown: blue contour: mostly lavas, pyroclastic density currents and block fans; green: July 1968 pyroclastic density current deposits and epiclastics (westward), and 1975 pyroclastic density current deposit on NW; Violet: epiclastic-rich valleys and fans. Below, two comparative topographic profiles are displayed in N–S (left) and W–E (right) directions (their location is indicated by yellow dash lines on the maps). The black arrows that point down indicate the intersection of both profiles. Instituto Geográfico Nacional provided the Arenal digital elevation model based on 2014–2019 images. Local coordinates EPSG 5456 Lambert North

Full size image

As it is given in Table 1, 70% of the total dense rock equivalent volume corresponds to the lava field area, where most of the new erupted material accumulated. The other 30% was accumulated mostly on the ample fan farther the lava field to the west (Fig. 6), where most of the pyroclastic density currents from 1968 were deposited (isopachs in Alvarado et al. 2006), and a thick pile of reworked material from the lavas was deposited as well. The Tabacón valley on the north-western flank was filled by pyroclastic density currents from 1968, 1975 and 1993, adding some 0.5% of the volume erupted (Fig. 6). The erosion of the newly deposited lava flows, especially triggered by heavy rains and flash floods, deposited thick narrow fans along the Cedeño lake basin (north), Agua Caliente river and tributaries (southwest), Guillermina creek (northeast), Calle de Arena creek and its tributaries and Chato creek (east) (Fig. 6), and many other smaller creeks all around the volcano (Zapata and Soto 1991). Part of this reworked material is evidently not from the latest eruption, coming down the slopes by erosion caused by bombardment, deterioration of the upper cone vegetation and a renewal of the erosive pattern. It is very difficult to assess the volume eroded and redeposited, as well as to calculate the apparent volume added by an increase in the porosity of the new deposits compared to the older ones eroded. If we take into account that 35% of the total new volume of the cone is from pyroclasts (~ 50% of the total) and epiclastics (~ 50% of the total), and from those epiclastics less than one-third is older material (cf. Zapata and Soto 1991) with an increase of 50% of its original volume, we have no more than 3% of the total volume as a bulking up of the older material reworked and redeposited as apparently new. This is one-third of the uncertainty and therefore probably included into.

Discharge rates

The evolution of the discharge rate is shown in Fig. 5a. The displacement of eruptive activity from crater A to C reduced the discharge rate by a factor of 2 (Wadge 1983). To explain this, Wadge (1983) modelled a plumbing system in which craters A and C were located above extensive eastward-dipping conduits that converge towards a point vertically below the crater D and attributed the drop in the discharge rate to the difference in static pressure generated by the additional weight of the 400 m of magma between the craters. Such inverse correlation between the mean effusion rates and the elevation of the effusive vents was observed at Stromboli (Ripepe et al. 2017) and is driven by the magma column stored in the central conduit system above the lateral vent (Ripepe et al. 2015, 2017). It means that a steady supply of magma from depth is needed to maintain the long-lasting eruption at Arenal in addition to the gravity-driven discharge dynamics (Ripepe et al. 2015), as supported by petrological and geochemical data (Streck et al. 2002; Gill et al. 2006) and SO2 estimates (Williams-Jones et al. 2001; Wade et al. 2006).

Petrology and rheology

The erupted products of Arenal range from basaltic to dacitic in compositions. While the SiO2 content for tephra ranges from ~ 49 to 63 wt % that for the lava spans within a narrower range of ~ 54 of 56 wt % (Borgia et al. 1988; Bolge et al. 2006). There is consensus over a mid-crustal (~ 10 km deep) provenance of Arenal magmas (e.g. Cigolini and Kudo 1987; Cigolini 1998; Streck et al. 2005; Ryder et al. 2006). This is also consistent with geodetic observations (Wadge 1983; Ebmeier et al. 2010, 2014; Mora et al. 2013) and volatiles in melt inclusions that show degassing coupled with magma fractionation and ascent from ~ 2 kbar (~ 7 km; Wade et al. 2006). This magma source is also confirmed by seismic tomography in which a negative velocity anomaly beneath Arenal has been observed at a depth of between 4 and 9 km (Núñez et al. 2020).

Phenocryst-rich basaltic–andesites erupted during the 1968–2010 eruptive phase have an SiO2 content in the range 54–57 wt % (Reagan et al. 1987; Streck et al. 2005; Ryder et al. 2006). Tephra deposits from the opening explosive phase (29–31 July 1968) have some of the highest SiO2 contents (56–57 wt %) and indicate slightly zoned reservoir at the beginning of the eruption (Alvarado et al. 2006). Despite this, small compositional range mineral compositions record a complex magmatic history during the 1968–2010 eruption (Parat et al. 2014). Multiple geochemical models have been proposed to explain the source and evolution of Arenal magmas (e.g. Reagan et al. 1987; Cigolini 1998; Streck et al. 2002, 2005; Ryder et al. 2006; Szramek et al. 2006). The models agree on a multi-stage evolution of the magmatic system starting from a stratified magma reservoir perturbed by a deeper basaltic intrusion, leading to the onset of the eruption in 1968. However, the models differ on the mechanisms of magma generation and evolution during the remainder of the eruption. Reagan et al. (1987) initially proposed a three-stage process of fractional crystallisation, magma mixing and fractional crystallisation plus recharge. Instead, Cigolini (1998) defined an evolution implying assimilation, fractional crystallisation, recharge and, finally, crystal fractionation. Streck et al. (2002, 2005) concluded that Arenal magmas resulted from multiple mixing stages. Finally, Ryder et al. (2006) proposed a simpler two-stage model, subdividing the second one into four stages, invoking only fractional crystallisation and recharge to explain most of the observed variations in bulk composition. Regardless of the number of stages proposed, the first half of the 1970s marked a change in the system dynamics, with the ascent rate (Szramek et al. 2006) and discharge rate (Wadge 1983; Wadge et al. 2006) both decreasing. In parallel, in 1974 the eruptive activity shifted from crater A to crater C (Wadge 1983; Fig. 1b). Ryder et al. (2006) proposed that early 1971 marked the shift from closed to open-vent activity.

Different works studied the long-lasting and steady nature of the 1968–2010 eruption (e.g. Wadge 1983; William-Jones et al. 2001; Streck et al. 2002), even though none of them has covered the final stage (from 2000 to October 2010). Streck et al. (2002) inferred that Arenal is underlain by a continuously active, small-volume magmatic reservoir. They argued that the reservoir was maintained in a quasi-steady state by basalt recharge over several decades. This hypothesis is supported by the magma budget estimates carried out by Wadge (1983) and Wadge et al. (2006), SO2 fluxes obtained by William-Jones et al. (2001) and the geodetic studies of Ebmeier et al. (2010) and Mora et al. (2013). Gill et al. (2006) pointed out that the period from 1968 to 1970 involved the eruption of increasingly mafic magma. However, from 1970 to 2000, the steady phase was dominated by recharge by a chemically different and more mafic magma that probably initiated and sustained the eruption up to its end (Gill et al. 2006).

Kozono and Koyaguchi (2009) obtained a viscosity for a parental hot (950–1050 °C) and hydrous (4 wt. % H2O) magma (crystal-liquid mix) in the range of 102.5–104 Pa∙s. To calculate this, they used a combined model from Hui and Zhang (2007) and Costa (2005) along with geochemical and petrological data from Reagan et al. (1987), Streck and Wacaster (2006) and Wade et al. (2006). In contrast, Bertolino et al. (2009) obtained viscosities in the order of 106 Pa∙s for a melt with crystallinity 42–54 vol. % at or near the surface and under dry conditions by using a modified Tamman–Vogel–Fulcher (M-TVF) equation combined with the modified Einstein–Roscoe (MER) equations. The viscosities were explained by Wade et al. (2006, Fig. 2), who showed magma differentiation from a basaltic magma (4 wt. % H2O) at 0.2 GPa (~ 5 km) to a basaltic andesite to andesitic magma (1 wt. % H2O) at < 0.02 GPa (< 1 km). Crystallinity at or near the surface could also reach ~ 55% (Reagan et al. 1987; Szramek et al. 2006), also contributing to the increase in viscosity of the mixture (Cigolini et al. 1984; Bertolino et al. 2009).

What we have learned from a geophysical perspective?

Seismicity at Arenal: 1968–2010

Seismological investigations at Arenal began shortly after the beginning of the 1968 eruption, when Matumoto (1968) deployed instruments from 6 to 14 August 1968, at La Fortuna and La Palma, 6.8 km east and 3.6 km north–north-east from the active crater, respectively. This installation resulted in counting of microearthquakes (typically magnitude 2–3 and no greater than 3.5) and made the first description of microtremors. Later, Minakami et al. (1969) installed the first two large aperture arrays (three stations each) at the same locations. The first array at La Fortuna had distances between seismometers of 1.4 km to 3.4 km and started recording on 25 and 29 August 1968, while the array at La Palma had distances between seismometers of ~ 1 km and began recording on 29 September 1968. During two weeks of observations, Minakami et al. (1969) described volcano-tectonic, long-period and tremor events. Minakami et al. (1969) also made the first description of an explosion quake, as recorded on 16 September 1968 at 21:49:34, which had an impulsive positive onset at all stations that they explained by a compressive source.

Numerous studies came over the following five decades, focusing on understanding Arenal's seismicity and its relation with the eruptive activity. From 1968 until the late 1990s, most of the studies described the waveform and the spectral content of the seismic signals without accounting for the temporal evolution of their characteristics (Minakami et al. 1969; Matumoto 1976; Alvarado and Barquero 1987; Morales et al. 1988; Melson 1989; Barboza and Melson 1990; Barquero et al. 1992; Métaxian et al. 1996; Mora 1999a; Soto et al. 1998). However, from 1998 to 2010, seismological studies evolved towards more detailed descriptions and analyse by introducing time–frequency methods, revealing a wide variety of patterns (Benoit and McNutt 1997; Alvarado et al. 1997; Hagerty et al. 1997, 2000; Mora 2003; Lesage et al. 2006).

Álvarez et al. (2009) and Cortés et al. (2009, 2021) used descriptive databases of Arenal seismicity to develop and test automatic detection, classification and recognition systems of volcano-seismic events. Such systems are now more efficient and prove to be important tools for monitoring, especially in the case of a seismic crisis when a large number of events must be processed in a short amount of time (e.g. Benítez et al. 2007; Beyreuther et al. 2008; Arámbula-Mendoza et al. 2018).

Classification of seismic events

The characteristic Arenal seismicity includes five types of event that have since been known to characterise such persistently active open-vent systems: (1) explosion quakes, (2) long-period events, (3) harmonic tremor, (4) spasmodic tremor and (5) volcano-tectonic events. This has resulted in Arenal becoming the benchmark for studies and classification of seismic events at other open-vent systems such as Karymsky, Sangay, Fuego, Colima and Merapi (e.g. Johnson et al. 1998; Lees and Ruiz 2008; Lyons et al. 2010; Arámbula-Mendoza et al. 2011; Budi-Santoso et al. 2013). Here, we define each of these signals as witnessed in the Arenal record.

Explosion quakes

These signals are low frequency and concomitant with the explosions (Fig. 7a and b). An air-shock wave was frequently recognised several seconds after the P-wave onset, which was a characteristic feature from 1984 up to the late 1990s (e.g. Barquero et al. 1992; Métaxian et al. 1996; Alvarado et al. 1997; Hagerty et al. 2000). As also observed since at, for example, Pavlov, Stromboli, Sakurajima and Popocatépetl (e.g. Garcés et al. 2000; Ripepe et al. 2001a; Morrissey et al. 2008; Mendo-Pérez et al. 2021). The explosions were recorded up to a very long distance of 92 km (Barquero 1988), and the air-shock wave caused window vibrations and an audible boom up to several kilometres, even at Observatorio Sismológico de Arenal y Miravalles headquarters located 23 km west of the volcano. Métaxian et al. (1996) and Hagerty et al. (1997, 2000) observed that every explosion was recorded as an identical waveform, which is also characteristic at other open-vent volcanoes such as Karymsky, Erebus and Yasur (e.g. Johnson et al. 1998; Rowe et al. 2000; Battaglia et al. 2016). From 1998 up to 2010, this typical audible boom progressively disappeared as the explosive activity waned.

Fig. 7
figure 7

Examples of typical seismic events of Arenal. The left column depicts the normalised velocity waveforms and the right column the corresponding spectrograms. a Explosive event with air-shock wave. b Explosive event without air-shock wave. c “Whoosh” event. d Harmonic tremor. e Spasmodic tremor

Full size image

Alvarado and Barquero (1987), Morales et al. (1988) and Barquero et al. (1992) recognised that long-period events and explosion quakes could have harmonic tremor-like codas (Fig. 7c). These complex signals correspond to small explosions (so called “whooshes”) followed by rhythmic gas emissions accompanied by harmonic tremor (so called “chugs”) (Melson 1989; Barboza and Melson 1990; Benoit and McNutt 1997). The long tremor-like coda of the chugging used to have a fundamental frequency gliding from about 2 to 3 Hz (Benoit and McNutt 1997; Lesage et al. 2006).

Long-period events

Long-period events are discrete transients with an emergent onset and energy distributed at frequencies in the range 0.5–5 Hz (e.g. Chouet 1996; Lyons et al. 2014; McNutt and Roman 2015) and have also since being recognised at other open-vent systems such as Stromboli, Soufriere Hills (Montserrat) and Etna (e.g. Acernese et al. 2004; Neuberg et al. 2006; Cannata et al. 2015). These types of signals were little described in the specific case of Arenal. However, studies agree that they are very similar in frequency range (1–3 Hz) and waveform to explosions, except that their records do not show the acoustic phase. Consequently, there are no fundamental differences in source mechanisms for long-period events and explosion quakes (Métaxian et al. 1996; Hagerty et al. 2000; Lesage et al. 2006).

Harmonic tremor

Harmonic tremor is a persistent ground vibration lasting from several minutes to several hours and even days. It has a regular wave appearance in the time domain and exhibits a set of narrow, evenly spaced frequency peaks with a fundamental frequency and its harmonic overtones (e.g. Konstantinou and Schlindwein 2003; McNutt and Roman 2015; Saccorotti and Lokmer 2021). It is an extremely common signal at open-vent systems such as Semeru, Sakurajima, Popocatépetl and Shinmoedake (e.g. Schlindwein et al. 1995; Maryanto et al. 2008; Arámbula-Mendoza et al. 2016; Natsume et al 2019) and has been the subject of a wide range of studies to model its source processes (e.g. Julian 1994; Chouet 1988; Julian 2000; Hellweg 2000; Bercovici et al. 2013; Dmitrieva et al. 2013; Montegrossi et al. 2019; Takeo 2021).

Harmonic tremor was the most conspicuous signal at Arenal and the most studied. Minakami et al. (1969), Alvarado and Barquero (1987), Alvarado et al. (1988b), Barquero et al. (1992) and Métaxian et al. (1996) reported fundamental frequencies as low as 0.5 and 0.65 Hz. Studies conducted after 1997 do not report such low frequencies. Instead, they are fairly consistent in reporting fundamental frequencies at ~ 0.9 Hz (Hagerty et al. 1997, 2000; Mora 1999a, 2003; Lesage et al. 2006) or in the range of 1.8–1.9 Hz (Alvarado et al. 1997; Benoit and McNutt 1997; Hagerty et al. 1997, 2000; Soto et al. 1998; Mora 2003; Lesage et al. 2006). Furthermore, studies consistently reported a peak at 1.8–2 Hz (which is an even overtone of the 0.9 Hz frequency) as the largest, regardless of the type of instrumentation or observation site (Fig. 7d).

Spasmodic tremor

Spasmodic tremor is an irregular, sometimes pulsating signal, characterised by its long duration (hours to days) and broad spectrum reaching frequencies above 5 Hz (e.g. Konstantinou and Schlindwein 2003; McNutt and Roman 2015; Saccorotti and Lokmer 2021). This type of tremor has been observed, for example, at Krafla, Galeras, Deception Island, Popocatépetl, Kelud (e.g. Brandsdóttir and Einarsson 1992; Gil-Cruz 1999; Almendros et al. 1997; Arámbula-Mendoza et al. 2016; Hidayati et al. 2019) and has been associated with magma intrusion, long-period event swarms, sustained water steam phases, ash venting and volcano-tectonic events bursts. At Arenal, this non-harmonic tremor is usually in the range of 1–6 Hz, with the record being dominated by multiple uncorrelated frequency peaks (Alvarado et al. 1997; Mora 2003; Lesage et al. 2006) (Fig. 7e).

Volcano-tectonic events

Volcano-tectonic events are high (2–20 Hz) frequency events with distinguishable P and S phases, generated by brittle fracture and movement along faults (e.g. Roman and Cashman 2006; White and McCausland 2016; Saccorotti and Lokmer 2021). According to Roman and Cashman (2006), volcano-tectonic events have been recorded in a variety of volcanic settings and types of magma involved in eruptions, such as: Mount Unzen eruption in 1991 (extensional graben/arc, dacite), Soufriere Hills eruption in 1995 (compressive/arc, andesite), Miyake-jima eruption in 2000 (compressive/triple junction, basalt–basaltic andesite) and Etna eruption in 1991 (mixed compression/extension, basalt).

Little is known about volcano-tectonic seismicity at Arenal. A few studies have described volcano-tectonic swarms prior to the eruptions of 1968 and 1975 (Matumoto 1976; Barquero et al. 1992; Alvarado et al. 2006), as well as prior to some pyroclastic density currents (Barquero et al. 1992; Alvarado et al. 1998; Arroyo et al. 1999, 2000; Alvarado and Soto 2002). Hagerty (1998) observed that much of the seismicity reported as tectonic was indistinguishable from the explosion quakes, leading him to think that it was possible that much of the reported seismicity was due to explosions. Furthermore, Hagerty (1998) also mentioned that volcano-tectonic activity was rare in Arenal, which is also reinforced by the fact that there are practically no reports of it in the Arenal literature. Alternatively, it is possible that this seismicity was of low magnitude and difficult to recognise due to the great amplitude of low-frequency seismicity caused by the eruptive activity.

We obtained the volcano-tectonic event catalogue from the Observatorio Sismológico y Vulcanológico de Arenal y Miravalles, which extends from 2001 to 2021. Events within a range of 3 km from the summit crater and with a magnitude greater than 1.5 were selected. The events with large uncertainties were also filtered out. The remaining locations are shown in Fig. 8. We can observe that most seismic activity is located on the western and southern flanks of the volcano at shallow depths (< 6 km) and with magnitudes less than 3.5. This seismicity could correspond to: 1) activity of regional faults delimiting the graben on which Arenal is built, as well as other faults crossing the cone, 2) a magma feeding system that may be slightly shifted towards the west, and/or a combination of both. A more detailed study of this seismicity would be required to assess these hypotheses for the origin of such volcano-tectonic events.

Fig. 8
figure 8

Volcano-tectonic seismicity at Arenal from April 2001 to August 2021 (M ≥ 1.5). The colour code indicates the date (year-month) of the events and the size of the symbols depends on their magnitude. Instituto Geográfico Nacional provided the Arenal digital elevation model based on 2014–2019 images. Local coordinates EPSG 5456 Lambert North

Full size image

Precursory seismicity

The precursory seismicity of the 29 July 1968 eruption consisted of volcano-tectonic swarms that were felt by the inhabitants (Table 2). The deployment of the first seismic network in 1974, allowed the instrumental observation of an eruption on 17 June 1975, and some precursor swarms, for the first time in Costa Rica (Matumoto 1976; Barquero et al. 1992). However, it was not until 1998, when Alvarado et al. (1998) and later Alvarado and Soto (2002) reported systematic precursory seismicity prior to pyroclastic density current events derived from crater wall collapse (Table 2). Barquero et al. (1992) and Soto et al. (1998) suggested possible correlations of tremor frequency variations to pre-eruptive phases or increasing explosivity. Table 2 summarises observations of precursory seismic activity before the main eruptive episodes.

Table 2 Summary of the precursory seismicity observed at Arenal at different stages of its activity
Full size table

Studies of seismic sources

Behaviour of volcanic tremor

The application of time–frequency processing to the analysis of harmonic tremor at Arenal revealed the phenomenon of gliding, in which the set of spectral peaks shifts as a function of time while maintaining their regular spacing (Benoit and McNutt 1997; Hagerty et al. 1997, 2000; Mora 2003; Lesage et al. 2006). Indeed, this phenomenon produced large variations (> 50%) of the fundamental frequency in minutes or tens of minutes (Benoit and McNutt 1997; Hagerty et al. 2000), which explains the different estimates of the fundamental frequencies on studies carried out in 1968–1996. Mora (2003) and Lesage et al. (2006) applied high-resolution, time–frequency methods that highlight other types of behaviour of harmonic tremor. Some of them are:

1) frequency jumps with either positive or negative increments.

2) tremor episodes with two simultaneous systems of spectral peaks subject to independent frequency gliding, denoting that two oscillation states may coexist for short periods, in which both even and odd overtones are observed.

3) progressive transitions from spasmodic to harmonic tremor.

4) absence of a clear and systematic relationship between the occurrence of explosions and tremors.

These complex features of Arenal’s tremor were not common at other volcanoes. However, their detailed study helped to understand the physical processes that produce this type of seismic vibration and was the basis for developing a conceptual source model for Arenal (Lesage et al. 2006; Valade et al. 2012; Almendros et al. 2014). This model will be presented in the discussion.

Location and mechanism of seismic sources

Minakami et al. (1969) derived the back-azimuth of an explosion quake recorded on 16 September 1968 at 21:49:34 using two tri-party large aperture arrays. The back-azimuth was consistent with the direction of the lowest crater (crater A, Fig. 1b). Moreover, Minakami et al. (1969) estimated an apparent P-wave velocity of 4.2 km s−1.

From the 1970s to 2000, the location of seismo-volcanic sources at Arenal was focused on explosive events by using clear P-wave arrivals combined with the acoustic air-shock wave phase (Alvarado and Barquero 1987; Alvarado et al. 1997; Hagerty 1998; Hagerty et al. 2000). After 2000, with the deployment of various experiments, more sophisticated location techniques were applied not only to explosive events but also to the tremor. The techniques were based on seismic arrays (Métaxian et al. 2002, 2009; Almendros et al. 2014), the inversion of the seismic moment tensor (Davi et al. 2010, 2012; Davi 2011) and coda wave interferometry (Snieder and Hagerty 2004).

Métaxian et al. (2002) developed a method based on a network of several triangular seismic arrays deployed around the volcano. They estimated the probability density function of the slowness vectors at each site, and by combining the probability density function of the back-azimuths, they obtained the probability density function of the source location. The position of the source was then defined as the maximum likelihood of this probability density function. They located the source of many tremor events, long-period events and explosions with this approach, and they defined in the horizontal plane, a seismogenic zone centred on the crater C (Fig. 9). Métaxian et al. (2009) calculated synthetic seismograms to quantify the scattering effect due to the complex topography and the shallow velocity structure. The results can be used to select the sites that minimise the topographic effects and to improve the precision of source location.

Fig. 9
figure 9

source by seismic triangulation (modified from Métaxian et al, 2002). a Seismogram of the located tremor episode. b Light brown triangles denote the positions of the triangular arrays. The probability density functions of the back-azimuth of propagation are represented in a polar diagram with 1° increments. The corresponding scale (in %) is displayed on the upper right corner of the map. The probability density function of the source position is displayed in colour with the scale in the lower right corner of the map. Its maximum gives the estimated source position (white cross)

Location of a tremor

Full size image

Almendros et al. (2012, 2014) carried out a detailed analysis of the wavefield recorded by a temporary array of 19 short-period seismometers with a 210-m aperture. They detected strong and time-varying differences in the amplitude of harmonic tremor recorded by receivers separated by some tens of metres. This complex pattern was not observed for spasmodic tremor or other types of events. They interpreted these observations as a phenomenon of interference between two or more components of the wavefield, which appears to be the first and unique description of this phenomenon for seismic waves. The long-period events and explosions were characterised by stable back-azimuth that pointed towards the crater with a fluctuation of ± 10°. In contrast, the harmonic and spasmodic tremors had more complex slowness vectors with back-azimuths that varied in an interval of ± 40° around the direction of the crater. These results have implications on the sources of events, as discussed next.

Davi et al. (2010) retrieved the source mechanism of an explosion on 14 February 2005 through moment tensor inversion. They inverted the waveforms recorded by nine three-component stations and made synthetic tests to study possible spurious single forces in the source representation. The resulting mechanism was mainly isotropic (87%), which corresponds to an explosive event producing a volume change of 68 m3 (assuming a shear modulus, μ, of 10 GPa) at about 200 m beneath the crater. Moreover, Davi et al. (2012) addressed, for the first time, the difficult issue of applying the moment tensor inversion approach to study the source mechanism of harmonic tremor. They obtained a constant mechanism for the entire duration (100 to 600 s) of all the 15 studied tremor events. The source was represented by a crack dipping 20° towards the north–north-east at shallow depth (~ 100 m).

Snieder and Hagerty (2004) also carried out one of the first applications of the coda wave interferometry method to volcanoes using data from Arenal. To do this, they deconvolved the harmonic tremor signal with air pressure records to obtain repetitive displacements of diffuse waves. They interpreted the slight changes in the deconvolved signals as due to a shift of about 15 m laterally in the source position. Since then, coda wave interferometry has been used to detect small velocity or structural changes in many volcanic structures (e.g. Brenguier et al. 2008; Obermann et al. 2013; Budi-Santoso and Lesage 2016).

Velocity models and site effects

A few studies have focused on Arenal's velocity models as well as topographic, path and site effects on the volcanic wavefield. Métaxian et al. (1996) used the arrival time differences between P and acoustic waves, as observed at four stations, to estimate a lower bound of the P-wave velocity as Vp > 1.3 km s−1. Hagerty et al. (2000) analysed records from a small linear array and obtained an estimation of the apparent P-wave phase velocity in the range of 2.5–4.6 km s−1. Leandro and Alvarado (1999) carried out seismic refraction profiling and vertical geoelectric surveys on the eastern and western flanks of the volcano (Fig. 10a and b). They obtained detailed 2D velocity models for depths of up to 300 m that revealed strong lateral and vertical velocity variations as well as marked changes in the layer thicknesses (Mora et al. 2006). They also found evidence of faulting, associated with low-velocity anomalies. Mora (2003) and Mora et al. (2006) estimated 1D velocity models at a depth of a few hundred metres using the spatial correlation method (SPAC) (Aki 1957) and data from two semicircular arrays located at 3.5 km NE (NESC) and 2.5 W (WSC) of crater C (Fig. 10a inset). The resulting simplified velocity model is summarised in Table 3. The velocities are consistent, down to a depth of 400 m, with those obtained by the seismic refraction models of Leandro and Alvarado (1999). Badilla and Taylor (2019) obtained a 7-km-deep magnetotelluric model parallel to the seismic refraction profile on the western flank of Arenal. They observed a resistivity range from 25 to 250 Ωm for the first 150 m that is consistent with the resistivity range of 25–1000 Ωm observed by Leandro and Alvarado (1999). Moreover, Badilla and Taylor (2019) also observed a major resistivity contrast at a depth of ~ 500 m, which is in good agreement with the velocity change at ~ 400 m obtained by Mora et al. (2006). Lesage et al. (2018) subsequently used the velocity models obtained at Arenal, together with those of other volcanoes, i.e. Masaya, Stromboli and Unzen (Métaxian et al. 1997; La Rocca et al. 2000; Sakuma et al. 2008), to define a generic velocity model to be used at volcanoes when none is available.

Fig. 10
figure 10

a Map showing the location of seismic stations deployed on the western flank of Arenal volcano during several experiments (coloured circles) and the track of the seismic refraction survey (brown line). The range of resonance frequency of each site, estimated by the H/V spectral ratio, is indicated by the colour of the circles. Grey contours lines every 10 m. Instituto Geográfico Nacional provided the Arenal digital elevation model based on 2014–2019 images. Local coordinates EPSG 5456 Lambert North. b Simplified 2D velocity structure obtained by seismic refraction survey (Leandro and Alvarado 1998). The stations of a linear array deployed in 1997 (Mora et al. 2006) are numbered in both panels

Full size image
Table 3 Simplified velocity model obtained with the Spatial autocorrelation method
Full size table

The strong variations in spectral peak amplitudes observed along the seismic arrays deployed in 1997 were attributed to effects produced by the recording sites by Mora (1999b, 2000, 2003) and Mora et al. (2001, 2006). This interpretation led Mora (2003) and Mora et al. (2006) to carry out detailed studies of these effects by combining the horizontal-to-vertical (H/V) spectral ratio method (Nakamura 1989) and S-wave theoretical transfer functions obtained from the seismic refraction surveys. They demonstrated that the amplification of wave amplitudes was related to the resonance of a surficial structure comprising pyroclastic and epiclastic deposits (with Vp < 0.9 km s−1) and occurs only when the impedance contrast with underlying layers with higher Vp is strong enough.

To support our review, we thus calculated the horizontal-to-vertical (H/V) spectral ratios for 48 sites located along the western flank of Arenal where previous seismic deployments have been conducted. We thus reused data from Hagerty et al. (2000, station WARN), Métaxian et al. (2002), Mora et al. (2001, 2006), Guerrier (2012), Almendros et al. (2012, 2014) and the reference station of RSN, ARE1 (Fig. 10a). We extracted segments (8 to 15 min) of background noise and followed the procedure of Mora et al. (2001, 2006) to calculate the H/V ratios. At each site, we obtained the resonance frequency of the shallow structure (Fig. 10a). We can observe that resonance frequency ranges between 1 and 2 Hz for most of the sites and that the highest frequencies were obtained towards the west of the area, where the very shallow layers become thinner (Fig. 10b).

Acoustic

Strong explosions produce acoustic waves that propagate in the atmosphere and can be transmitted to the ground and thus also detected by seismometers, such as at Karymsky (Russia) and Sangay (Ecuador) (Johnson and Lees 2000), Arenal (Hagerty et al. 2000), Stromboli (Ripepe et al. 2002) and Erebus (Johnson and Aster 2005) (Johnson et al. 2003; Arrowsmith et al. 2010). The audible sounds that accompanied the explosions at Arenal, such as the “whooshes” and “chugs”, were recorded by both permanent and temporary seismic and acoustic deployments. These data can be used, not only for locating sources (e.g. Ripepe and Marchetti 2002; Johnson 2004; Johnson et al. 2011), but also for understanding source processes, such as for the temporal evolution of either the flow front of a lahar (Johnson and Palma 2015) or the level of the lava lake (Goto and Johnson 2011; Richardson et al. 2014; Johnson et al. 2018). Acoustic arrays can also be used to characterise the type of acoustic sources for strombolian explosions, such as at Erebus (Johnson et al. 2008; Gerst et al. 2013) and Yasur (Iezzi et al. 2019; Fee et al. 2021) or for vulcanian explosions such as at Sakurajima (Johnson and Miller 2014; Kim et al. 2015; Fee et al. 2017) or Fuego (De Angelis et al. 2019), leading to a robust estimate of gas volumes and gas fluxes. Furthermore, acoustic arrays can provide an estimate of the crater geometry for extended sources, such as for the Santiaguito’s dome (Johnson et al. 2011) or for the vent at Etna (Watson et al. 2020).

In 1997, a linear array of seismometers and microphones was deployed at about 2 km from Arenal’s active crater (Garcés et al. 1998). These seismic and acoustic records displayed great similarity. In particular, during the episodes of harmonic tremor, seismic and acoustic spectrograms displayed evenly spaced spectral peaks (at the same frequencies) and identical frequency gliding (Garcés et al. 1998; Mora et al. 2009). The acoustic pressure associated with the explosions had amplitudes of up to 125 Pa at 2 km from the crater (Hagerty et al. 2000). However, the energy partitioning between seismic and acoustic signals was strongly variable for the explosions, although the delay between arrival times of the two phases was stationary. This is consistent with a stable source depth as observed at volcanoes such as Erebus, Karymsky and Sangay (Johnson et al. 2003), at Stromboli (Ripepe et al. 2001a) and Etna (Ripepe et al. 2001b). By analysing this delay and by modelling the acoustic waveforms, the depth of explosion sources at Arenal was estimated in the range of 10–100 m below the crater surface (Hagerty et al. 2000), which is consistent with the depth estimation of Alvarado and Barquero (1987) and Alvarado et al. (1997) which was less than 250 m beneath the crater.

We present the results obtained from broadband acoustic measurements carried out from 11 to 22 February 2005 with the aim of assessing the type of activity and the gas volume expelled at Arenal volcano during this period. The acoustic and seismic instrumentation was deployed at 1.7 km from crater C (Fig. 1b) with a direct view to the vent. It comprised a microphone (Bruel Kjaer model 4193, 0.1–5 kHz frequency range) and its pre-amplifier (Bruel Kjaer model 2669), a vertical 1 Hz seismometer (L4C-1C) and a datalogger (REFTEK 130) with a sampling frequency of 1000 Hz, restricting the frequency band to 0.1–500 Hz. During the 11 days of observation, we only recorded five small explosions but numerous small pressure transients at a mean rate of seven events every six hours (Fig. 11a). Manual extraction was necessary due to the low energy of the small pressure transients, but allowed us to define these two seismo-acoustic event types at Arenal and compare them with similar waveforms recorded elsewhere.

Fig. 11
figure 11

a Number of small pressure transients per 6-h window, starting on 11 February 2005 at 18:40:16.97 UTC. Explosions are indicated by a red triangle. b Stacked waveform of small pressure transients recorded within the window starting at 12 February 2005 at 06:40:16.97 UTC (blue) and its corresponding synthetic acoustic pressure based on a bubble volume mode of a large bubble prior to its breaking and using the assumption of the magma layer thickness above the bubble to be 8 cm at equilibrium and the conduit radius equal to 2.5 m (Vergniolle and Brandeis 1996). c On the left is the ground velocity of one of the strongest explosions and its spectrogram (right). d Acoustic pressure (Pa) of the same explosion (left) and its spectrogram (right). The inset shows a zoom of the acoustic pressure along the window indicated by red dashed lines. Both waveforms were low-pass filtered at 10 Hz and the spectrograms estimated using a 10-s window and 50% overlap. Note the well-marked harmonic peaks following the explosion, which qualifies this explosion as a “whoosh”

Full size image

February 2005 explosions

The five recorded explosions, whose maximum radiated acoustic pressures reached 7.6 ± 4.3 Pa, were clearly detected by the seismometers (Fig. 11c and d). The strongest explosion was recorded on 14 February 2005 and generated an ash plume lower than 1 km. The moment tensor of this explosion was also studied by Davi et al. (2010).

The acoustic waveforms and the intensity of the explosions are similar to those recorded at Erebus and Yasur volcanoes (Johnson et al. 2008; Iezzi et al. 2019; Fee et al. 2021), suggesting that the source of the sound at Arenal is also produced by a monopole, i.e. by a sudden variation of the mass flux. The Arenal's waveforms are also similar to those recorded at volcanoes with intermediate to evolved magma composition, such as Karymsky and Sangay (Johnson 2003), Fuego (Lyons et al., 2010; Díaz-Moreno et al. 2020), Sakurajima (Johnson and Miller 2014; Kim et al. 2015) and Santiaguito (Sahetapy-Engel et al. 2008; De Angelis et al. 2016; Carter et al. 2020). However, they differ from those at Stromboli (Vergniolle and Brandeis 1996; Ripepe and Marchetti 2002; Ripepe et al. 2002), Etna (Vergniolle and Ripepe 2008; Spina et al. 2015), Erta Ale (Bouche et al. 2010) and Villarica (Ripepe et al. 2010; Johnson et al. 2018), which are either symmetrical or asymmetrical in favour of the rarefaction. They also differ in their maximum acoustic pressure larger at Arenal by at least a factor 2 than at these basaltic volcanoes such as Stromboli, Etna or Erta Ale. However, the maximum radiated acoustic pressure is similar to strombolian explosions of intermediate strength at Yasur and Erebus (Johnson et al. 2008; Spina et al. 2016; Witsil and Johnson 2018; Iezzi et al. 2019; Fee et al. 2021) and often much smaller than at the volcanoes with intermediate to evolved magma composition, such as Karymsky, Sakurajima and Santiaguito. The duration of the explosions also varies, from very short, less than a few seconds for Erta Ale (Bouche et al. 2010) or some gas-rich explosions at Stromboli (Vergniolle and Brandeis 1996; Ripepe and Marchetti 2002; Ripepe et al. 2002) or Etna (Vergniolle and Ripepe 2008; Spina et al. 2015) to duration exceeding 10 s for ash-rich explosions at Stromboli (Ripepe and Marchetti 2002) or at Yasur (Spina et al. 2016) and at volcanoes with more evolved magma composition, such as for the volcanoes mentioned above (Marchetti et al. 2009). The largest Arenal’s explosion of our period lasts probably between 10 and 15 s, but we cannot discuss this parameter or the type of acoustic source after the main pulse due to both the use of a sole infrasonic sensor and the low signal-to-noise ratio.

The robust assumption of a monopole source, valid for any type of explosion as long as the sound propagation is linear (Woulff and McGetchin 1976; Pierce 1981; De Angelis et al. 2019), is used to estimate the gas volume and gas flux during the main pulses of the infrasonic waveforms. The maximum gas flux and gas volume are obtained by one or two successive integrations of the acoustic pressure (Johnson et al. 2008; Johnson and Miller 2014; Iezzi et al. 2019; Fee et al. 2021). By applying this approach, we find that the gas volumes for the five explosions recorded here are between 1900 and 6000 m3 at atmospheric pressure. This is consistent with gas volumes estimated at other sites from infrasonic records, 5000–80 000 m3 at Yasur (Iezzi et al. 2019) and 1000–24 000 m3 at Erebus (Johnson et al. 2008). Another infrasonic modelling, based on the bubble volume mode prior to bursting, had been developed to estimate the gas volume for strombolian explosions (Vergniolle and Brandeis, 1996) (see details below). Values of gas volumes have been estimated at 36–700 m3 for Erta Ale (Bouche et al. 2010), 20–35 m3 (Ripepe and Marchetti 2002) and 2–100 m3 (Vergniolle and Brandeis 1996) for Stromboli, 750–7500 m3 for Etna (Vergniolle and Ripepe 2008) and 4700–30 000 m3 for Shishaldin (Vergniolle et al. 2004). The gas volumes produced by Arenal's explosions, 1900–6000 m3, are hence similar to the moderate to strong strombolian explosions elsewhere. Moreover, the maximum gas fluxes of the five Arenal’s explosions, 6200 and 20 000 m3 s−1, are also very similar to those estimated during strong explosions at Yasur, 3000–60 000 m3 s−1 (Iezzi et al., 2019), at Erebus 7000- 70 000 m3 s−1 (Johnson et al. 2008) and 4×104 m3 s−1 at Etna (Diaz-Moreno et al. 2019). The maximum gas fluxes on an hemispherical bubble head can also be estimated from the model based on the bubble volume mode prior to breaking by using the equilibrium radius of the bubble and its maximum radial velocity (Vergniolle and Brandeis 1996; Vergniolle and Ripepe 2008; Vergniolle et al. 2004), giving typical maximum gas fluxes of 510 m3 s−1 at Erta Ale (Bouche et al. 2010), 680 m3 s−1 at Stromboli and 1.2×105 m3 s−1 at Shishaldin. It is also encouraging to see that the maximum gas flux estimated at Etna by using the bubble volume mode, 2.6×104 m3 s−1, gives a very similar result to the single integration of acoustic pressure when assuming a monopole source, 4×104 m3 s−1 (Diaz-Moreno et al., 2019), showing that these two methods are robust. The maximum gas fluxes recorded during Arenal’s explosions, 3000–60 000 m3 s−1, are well within the range of values obtained during the strombolian explosions of Stromboli, Etna and Shishaldin. The vulcanian explosions at Sakurajima produce a maximum gas flux equal to 6–15×105 m3 s−1 (Kim et al. 2015; Fee et al. 2017), a value significantly larger than that of Arenal.

February 2005 small pressure transients

The small pressure transients, whose maximum radiated acoustic pressures were only 1.8 ± 0.80 Pa (Fig. 11b), were not accompanied by any seismic signal. Field observations during the experiment revealed an association between one small pressure transient with a brown-coloured cloud emission, suggesting a degassing event that occurred ~ 3 min before the strongest explosion (14 February 2005). These low-level pressure transients are evocative of a small degassing activity, called puffing, initially described at Stromboli, and shown to produce small and very frequent infrasonic pulses, very well correlated with the tremor (Ripepe et al. 1996; Harris and Ripepe 2007). This puffing had been associated with the mild degassing of non-spherical bubbles or to weakly overpressurised Taylor bubbles, i.e. large and cylindrical bubble sometimes called slugs (Ripepe and Gordeev 1999; Ripepe et al. 2002; Colò et al. 2010; Tamburello et al. 2012). Puffing has been also shown to exist, based on thermal records at Stromboli, Etna and Yasur (Gaudin et al. 2017a). However, the maximum acoustic pressure of the small pressure transients recorded at 1.7 km from Arenal's crater, 1.8 Pa, is equivalent to a pressure of 9 Pa at 340 m, taking the geometrical attenuation into account. This is significantly larger than the few Pa recorded at the same distance at Stromboli. This difference in radiated acoustic pressure is likely to be associated with the larger magma viscosity at Arenal than that at Stromboli. Nevertheless, this suggests that the small pressure transients at Arenal might share a similar origin that gas puffs at Stromboli. Puffing activity had also be reported at other strombolian volcanoes, such as Masaya (Branan et al. 2008), Villarica (Gurioli et al. 2008) and Yasur (Bani et al. 2013; Spina et al. 2016). The number of puffs per minute is much smaller at Arenal, 0.019 event, than at Yasur from 1–2 events (Bani et al. 2013) to 40 (Spina et al. 2016), Villarica with 9 events (Gurioli et al. 2008), Masaya with 6–8 events (Brannan et al. 2008) and Stromboli with 20–136 events (Gaudin et al. 2017b).

The records of the small pressure transients were stacked (arithmetic average) over 6-h windows to increase the signal-to-noise ratios. The acoustic waveforms obtained after stacking show that, in contrast to explosions, these transients have a stronger rarefaction peak than its previous compressive peak (Fig. 11b). To our knowledge, this feature has rarely been reported at volcanoes with intermediate to evolved magma composition, excepted sometimes at Santiaguito for isolated explosions (Gottschämmer et al. 2021) and for pairs of explosions (Carter et al. 2020). These asymmetrical infrasonic waveforms are very characteristic of strombolian explosions, such as those produced at Stromboli, Shishaldin, Etna and Erta Ale volcanoes (e.g. Vergniolle and Brandeis 1994, 1996; Vergniolle et al. 2004; Vergniolle and Ripepe 2008; Bouche et al. 2010). However, acoustic waveforms of strombolian explosions are not always asymmetrical in favour of the rarefaction peak, such as at Erebus (Johnson et al. 2008; Witsil and Johnson 2018) or at Yasur (Iezzi et al. 2019). At volcanoes presenting asymmetrical acoustic waveforms towards a larger rarefaction peak than the compressive one, it has been proposed that the source of the sound is driven mainly by the oscillations of a large bubble prior to its breaking (Vergniolle and Brandeis 1994, 1996), a model validated by laboratory experiments for low-viscosity magma (Kobayashi et al. 2010). This model, which is appropriate for the lowest viscosities of Arenal’s magma (102.5–104 Pa∙s; Kozono and Koyaguchi 2009; Bertolino et al. 2009), providing that the acoustic waveforms are asymmetrical, is not valid for the highest range of viscosities (106 Pa∙s; Bertolino et al. 2009), as shown at Erebus (Gerst et al. 2013). For example, the bubble oscillations at Erebus are often shown to exist prior to explosions (half of the explosions) but as a small, but detectable, precursory signal, with a gas overpressure of 400 kPa on average, and not as the main source of the sound during the explosions, the gas overpressure being reduced at 100 kPa during the explosion (Gerst et al. 2013). But the magma viscosity at Erebus, 106.2 Pa·s (Le Losq et al. 2015), is markedly the largest among magma viscosities at other strombolian volcanoes. The exact value of the cut-off viscosity, below which the dominant source of infrasound is related to bubble oscillations and not to bubble bursting, is totally unknown and is a very complex problem to solve. Furthermore, if the source of the sound is a monopole for the entire range of viscosities (Kozono and Koyaguchi 2009; Bertolino et al. 2009), the integration of acoustic pressure leads to first an increasing gas flux, reaching its maximum when the acoustic pressure passes the zero between the compressive and the rarefaction peaks. After this point, the integration of a negative acoustic pressure leads to decrease the gas flux. The gas flux becomes zero at the end of the explosion, only when the infrasonic waveform is symmetrical or asymmetrical in favour of the compressive peak as for the Arenal’s explosions. However, the marked asymmetry of the acoustic waveforms of the small pressure transients towards rarefaction leads to a negative value of the gas flux before approaching the peak of the rarefaction wave. Reaching a negative gas flux prior to the end of the explosion would imply that the volcanic gas would return in the conduit after its expulsion, which is totally unrealistic. This does not favour the method of integration of acoustic waveforms to obtain gas volumes. The lack of an acoustic array does not allow the clear determination of whether the waveform asymmetry is a path or a source effect. However, we note that the acoustic waveforms for explosions do not present these characteristics, suggesting that the asymmetrical acoustic waveforms of the small pressure transients result mostly from the source. By only considering the first part of the acoustic waveforms corresponding to a positive gas flux, we obtain an estimate of the gas volume at atmospheric pressure of 4500 ± 3300 m3 per transient and 7700 ± 11 000 m3 s−1, for the maximum gas flux.

Alternatively, when using the bubble volume mode model (Vergniolle and Brandeis 1996), the fits between the synthetic and recorded waveforms are very good (Fig. 11b). This model, which mainly depends on gas volume, overpressure and magma viscosity, requires knowledge of the thickness of the magma layer above the oscillating bubble, a thickness which is assumed to be similar to the average size of ejecta. Unfortunately, no measurement of the scoria dimensions was available at the time of the acoustic recording. We thus choose here a value between that used for a strombolian bubble at Stromboli, 2 cm (Vergniolle and Brandeis 1996), 4–5 cm (Ripepe and Marchetti 2002) and 6–51 cm with a mean of 12 cm and a standard deviation of 5.4 (Bombrun et al. 2015), and those at Etna (10 cm, Vergniolle and Ripepe 2008), Shishaldin (15 cm, Vergniolle et al. 2004) and Erta Ale (22 cm, Bouche et al. 2010). While the choice of this parameter has an impact on the final results for gas volume and gas overpressure (Vergniolle et al. 2004), using a value of 8 cm gives the best match between synthetic and measured acoustic waveforms, for all the small pressure transients. In this case, the model gives a typical gas volume and gas overpressure in the conduit of 220 ± 160 m3 and 0.5–1.3×105 Pa, respectively, based on all the small pressure transients. The corresponding gas volume of the small pressure transients at atmospheric pressure is 460 ± 270 m3. Decreasing the film thickness by a factor 1.5 increases the gas volume in the conduit by a factor 1.6 while decreasing the gas overpressure by a factor 3 (Vergniolle et al. 2004). This is similar to the volume for passive gas puffing at Stromboli, 19–211 m3 with an average of 48 m3 (Table 6 in Harris and Ripepe 2007), 40–150 m3, i.e. equivalent to a gas mass of 10–30 kg (Ripepe et al. 2008), 1.4–21.5 m3 from an UV camera (Tamburello et al. 2012), 0.1–10 m3 from thermal camera (Gaudin et al. 2017b) and slightly less than at Masaya, 1000–10 000 m3 (Branan et al. 2008).

The gas volume estimates of the small pressure transients, 4500 ± 3300 m3 (for the integration method), are within the lowest range of the gas volume estimated for Arenal’s explosions, despite very different acoustic waveforms and maximum radiated acoustic pressure (Fig. 11). It is likely that the difference in radiated acoustic pressure between the two types of events, by a factor 4, results from a difference in gas overpressure, suggesting that the small pressure transients are due to passive degassing and hence are produced at a very shallow depth, such those associated with gas puffing at Stromboli (Ripepe and Gordeev 1999), while explosions are associated with active degassing and hence have a deeper origin. The lack of seismic signal associated with the small pressure transients, a feature shared by both the outgassing events at Santiaguito (Lamb et al. 2019) and for some mildly to moderately explosive strombolian explosions, such as at Erta Ale (Bouche et al. 2010), Stromboli (Vergniolle and Brandeis 1996) or Yasur (Kremers et al. 2013), reinforces this interpretation.

Gas volumes expelled during explosions had been estimated by various methods, such as using infrasound, open-path Fourier transform infra-red, radar, visible videos and thermal videos, or radiometers which measure the average temperature in a field of view (see details in Vergniolle and Gaudemer 2015). Estimates from infrasound, although often lower than those obtained from other techniques, are only considering the overpressurised gas flow, whose variations are able to radiate infrasound. A steady gas flow cannot be measured directly on infrasonic records (Johnson 2003) or on seismic records. While infrasonic records mainly measure the signature of the overpressurised gas, other techniques provide a global view of the total gas volume, including pressurised and non-pressurised gas (Vergniolle and Gaudemer 2015). Furthermore, the duration of the explosions is markedly smaller from infrasound, of the order of 1 s, than from thermal or open-path Fourier transform infra-red or radar records, > 10 s at Stromboli (Ripepe et al. 2002; Tamburello et al. 2012), Yasur (Kremers et al. 2013) and Erebus (Witsil and Johnson 2018), showing the existence of a significant steady gas flow at the end of these acoustic waveforms. The difference in the gas volume estimates, deduced from infrasound and from the height of the small lava fountains, had been explained at Etna by the existence of a very long bubbly wake trailing at the back of the bursting overpressurised Taylor bubble (Vergniolle and Ripepe 2008). The large magma viscosity of Erebus prevents such a mechanism of wake formation to occur in the conduit (Vergniolle and Métrich, in revision, this issue). In this case, the difference in the “small” gas volume estimated on infrasound and the inferred “large” size of cavity, left by the burst bubble, had been explained by an additional gas supply provided, either by a very large but not pressurised following bubble or a complex mixture of smaller bursting bubbles (Gerst et al. 2013). Therefore, we could conclude that the gas volumes deduced for Arenal’s explosions from infrasound are a low estimate, only accounting for the overpressurised gas, likely to have a deeper origin than the non-pressurised gas.

Doppler radar

Two field experiments were carried out with a Doppler radar in 2004 and 2005 (Donnadieu et al. 2005; Valade and Donnadieu 2011; Valade et al. 2012). Broadband seismometers were also deployed during the second campaign. Doppler radar VOLDORAD2 is a ground-based instrument, which measures the amount of solid ejecta and the distribution of their velocities as a function of time (Dubosclard et al. 1999, 2004). The records obtained during explosions display two simultaneous phenomena: a short-lived signal associated with the projection of ballistic blocks, and an ash plume of much larger mass with lower velocity and longer transit times. The joint analysis of seismic and radar data showed a large variability in both types of waveform and non-systematic relationships between them (Valade et al. 2012). Tephra emissions were observed during explosions and tremor, but also during seismic quiescence. No clear relationship was obtained between the parameters measured by the radar and the amplitude and frequency content of the coeval seismic signals (Valade et al. 2012).

Deformation

Several deformation studies have been carried out at Arenal, as summarised in Table S5 (Supplementary Material 4). Most of them agree that there is a long-term subsidence at Arenal (Sawdo and Simon 1969; Simon et al. 1970; Global Volcanism Program, 1979; Wadge 1983; Van der Laat 1988; Alvarado et al. 1988a, 2003; Soto 1991; Hagerty et al. 1997; Mora 2003; Mora et al. 2013). Some of the proposed explanations for the origin of the deformations observed at Arenal are:

  1. 1)

    the pressure variation in a shallow (< 4 km) magma chamber (Sawdo and Simon 1969; Simon et al. 1970; Global Volcanism Program 1979; Van der Laat 1988);

  2. 2)

    the response to the load produced by the lava field (Wadge 1983; Alvarado et al. 1988a, 2003, 2010; Soto 1991; Mora 2003; Müller et al. 2015);

  3. 3)

    surface processes in the volcanic edifice such as slow gravity-driven slide (Ebmeier et al. 2010) and frequent rock fall and shallow landslides (Ebmeier et al. 2014);

  4. 4)

    the combination of the load produced by the lava field and the compression caused to a weak, thick layer of Miocene sediments underlying a thin, stronger layer of Plio-Pleistocene lavas beneath Arenal (Müller et al. 2015); and

  5. 5)

    the consolidation of the weathered volcanic bedrock portion (weak and plastic) and incipient deformation due to basement spreading resulting from volcano loading (Alvarado et al. 2003, 2010).

Among these explanations, that of a shallow magma chamber is the least plausible as it is not supported by modelling (Mora et al. 2013) and the geological evidence instead points to a deep magma chamber (~ 10 km) (Wadge 1983; Alvarado et al. 1988b; Soto 1991; Mora et al. 2013) based on petrological and geochemical data (Reagan et al. 1987; Sachs and Alvarado 1996; Streck et al. 2005; Williams-Jones et al. 2001). The other sources are possible and even complementary, reflecting the complexity of the deformation occurring at Arenal. So far, no study has been able to determine the proportion of the contribution of each of them.

Discussion

Rose et al (2013) defined an open-vent system as the one with a continuous emission of magma-related products directly into the atmosphere, which takes place at a volcanic vent. Following this definition, Arenal became open-vent system in the early 1970s and remained that way until October 2010. The main features of the activity at this system were:

1) A lava pond that occupied the active crater from the second half of the 1980s to 1999 and evolved into formation of dome-like structures and occasional “hornitos” from 2000 to the end of the eruption.

2) Quasi-continuous effusion of lava flows.

3) Frequent explosive activity that continued at variable rates and that declined, progressively, during the last 10 years of the eruption.

4) Pyroclastic density currents by column collapse, as well as lava pond or lava front collapse.

5) Profuse degassing.

6) A low-frequency seismicity comprised by tremor (harmonic and spasmodic), explosion quakes and long-period events.

Persistent ascent of basaltic magma from the mantle generated this quasi-continuous activity, whereas complex processes of mixing, assimilation and differentiation explains the evolution towards basaltic–andesitic magma (e.g. Reagan et al. 1987; Cigolini 1998; Streck et al. 2005; Ryder et al. 2006; Wade et al. 2006). This quasi-continuous feeding is testified by a variable discharge rate that decayed, over four decades, from 0.93 m3 s−1 to 0.1 m3 s−1 (Wadge, 1983; Soto and Arias 1998; Wadge et al. 2006) and by the SO2 fluxes (William-Jones et al. 2001).

The lack of continuous long-term multi-parametric observations and the different temporal and spatial scales of the short-term geophysical experiments make it difficult to have a comprehensive record of the evolution of the open-vent eruptive phase in a geophysical sense. Nonetheless, multiple geophysical experiments and data obtained over decades of research allow us to address some key aspects about Arenal eruptive activity as we now develop in the discussion.

Mechanisms of explosive activity

There is no consensus on whether the Arenal explosive activity was purely strombolian or vulcanian. The typical strombolian regime consists of the ascent of large, buoyant gas bubbles bursting at the free surface (e.g. Wilson 1980; Vergniolle and Jaupart 1986; Vergniolle and Brandeis 1994, 1996; Taddeucci et al. 2015), associated with the lowest range of magma viscosities. One of the lowest estimates is that of Erta Ale, 22–27 Pa·s when using the model from Giordano et al. (2008) (Vergniolle and Métrich in revision, this issue), while the highest is probably that of Erebus, 106.3 Pa·s (Le Losq et al. 2015). The vulcanian regime, typically associated with viscous magma, results from the disruption of a conduit plug or dome that fails under a sufficiently high pressure gradient in the underlying magma (Clarke et al. 2015).

There are few descriptions of strombolian activity in Arenal. Barquero and Alvarado (1989a) and Dellino and Alvarado (1992) described juvenile vesicular pyroclasts typical of strombolian explosions obtained from eruptions that occurred in May 1988 and August 1992, respectively. However, none of them provide values for vesicularity.

Barquero (1980), Soto and Arias (1998) and Alvarado and Soto (2002) described infrequent formation of “hornitos” and intense cyclic (every ~ 40 s) glows that are more typical of a strombolian regime. In contrast, the textural analysis by Cole et al. (2005) of tephra fall produced by six explosions between 1987 and 1996 revealed a high proportion (50 to 71%) of blocky fragments compared to the rest that were partially fluid and partially blocky, fluid or vesicular. For Cole et al. (2005), the blocky fragments came from the destruction of the lava cap at the top of the conduit. This proportion of blocky fragments is comparable to the one of Sakurajima tephra (43 and 67%) reported by Miwa et al. (2013). Furthermore, Szramek et al. (2006) analysed ash samples obtained from three eruptions on 9 April 1990, 20 February 1996 and 1 September 1998, which they classified as strombolian, and obtained vesicularities of 4 vol. % for the first one and 10 vol. % for the other two. However, these values may be too low for strombolian eruptions fragments. At Stromboli and Villarrica, pyroclasts have vesicularities ranging from 43 to 73 vol. % (Lautze and Houghton, 2007) and 42–80 vol. % (Gurioli et al. 2008), respectively, while at Sakurajima vesicularities are < 10 vol. % and < 40 vol. % for blocky and vesicular fragments, respectively (Miwa et al. 2013). Therefore, it is plausible that the fragments analysed by Szramek et al. (2006) could have come from the fragmentation of the lava cap.

The textural analysis of the juvenile fragments of the 28 August 1993 pyroclastic density current carried out by Alvarado and Soto (2002) suggested that the conduit was vertically stratified, with an upper part formed by a colder and denser lava cap overlying a hotter magma with varying vesicularity. This view of the conduit reflects different degrees of cooling of a near-solidus lava pond prior to its fragmentation (Alvarado and Soto 2002). A similar texture-based stratification has been observed at Sakurajima which has a typical vulcanian regime during its open-vent eruption since 1955 (Miwa et al. 2013). Arenal and Sakurajima are also comparable in terms of viscosities. For hot and hydrous parental magma, viscosities range 102.5–104 Pa∙s and 103.2–104.5 Pa∙s (Kozono and Koyaguchi 2009), for Arenal and Sakurajima, respectively, whereas for the anhydrous case viscosities are 106 Pa∙s (Bertolino et al. 2009) at Arenal and ~ 107 Pa∙s (Miwa and Geshi 2012) at Sakurajima.

Melson (1989) and Barboza and Melson (1990) suggested that the more intense explosions were related to the formation of a lava cap, while smaller explosions (“whooshes”) with the extended “chugs” were related to a more open condition. Johnson et al. (2003) associated “whooshes” and “chugs” with a strombolian regime and suggested that they might be related to low yield-strength caps of rubble temporarily plugging vents. Barboza and Melson (1990) and Williams-Jones et al. (2001) reported an inverse correlation between tremor and explosive activity involving open (or partially open) and close conditions, respectively. Thus, Arenal was able to switch from vulcanian to strombolian regimes (“whooshes” and “chugs”) on a daily or hourly scale. This is consistent with Miwa et al. (2013), who observed, at Sakurajima, that the cap can be regenerated in 1–2 h. In addition, evidence of multiple venting in crater C could explain the coexistence of effusive and different types of explosive activity. Therefore, the explosive activity at Arenal may have been both strombolian and vulcanian (Johnson et al. 2003; Alvarado 2011) but in proportions that varied during all the explosive phases from 1984 to October 2010. Similar shift and/or coexistence from one regime to another has been observed at Sakurajima volcano during its open-vent activity period having a variable rate of vulcanian explosions but also fountains of lava at the Showa crater after chugging events, as well as strombolian eruptions and continuous ash emissions in Minamidake crater (Iguchi et al. 2020).

Understanding how and why eruptive regime changes within the same eruptive period or even from one eruption to the next at a given volcano is a fundamental but complex challenge (Cassidy et al. 2018). Different magma properties or parameters, such as the ascent and decompression rates, porosity, permeability, crystal and bubble content, rheology, viscosity and temperature, as well as the conduit and vent geometry control the eruptive explosivity and consequently the style (Cassidy et al. 2018). Bain et al. (2021), for example, proposed a model based on the decompression rates for Galeras volcano. High decompression rates (~ 0.167 MPa s−1), low-viscosity magma and SO2 fluxes (500–1000 tons day−1) lead to a rapid formation of low-permeability plug generating larger explosions at short repose times (tens of days), while low decompression rates (~ 0.0167 MPa s−1), higher viscosity magma and SO2 fluxes (7000–10 000 tons day−1) lead to a slow formation of a high permeability plug or dome generating small-volume explosions at long repose times (hundreds of days).

At Arenal, Szramek et al. (2006) estimated that the 1968 initial violent eruptive phase was associated with magma decompression rates between 0.0013 and 0.025 MPa s−1 and ascent rates of 0.05–0.9 m s−1. However, after the mid-1970s, the estimated decompression rate was < 0.0013 MPa s−1 and the ascent rate < 0.05 m s−1, causing Arenal to shift to a more effusive and less explosive activity (Szramek et al. 2006). At the same time, the discharge rate decreased from 2 m3 s−1 during 1968–1973 to 0.3 m3 s−1 thereafter (Wadge 1983; Wadge et al. 2006). Moreover, Szramek et al. (2006) proposed that a small change in magma ascent rate could generate differences between explosive and effusive activity, without changes in the magma composition. According to Bain et al. (2021), magma ascending at a slower rate has more time to transfer heat to the conduit walls than a magma ascending more rapidly, which may be expected to retain more heat when it reaches the shallow conduit. Hence, high temperatures contribute to lower viscosity in a rapidly decompressed andesitic magma and low temperatures contribute to increased viscosity in a slowly decompressed magma. This could be a plausible explanation for variations in explosive regimes and alternation between effusive and explosive periods at Arenal. Changes in the rheology of the upper part of the conduit and the lava pond, with frequent development of caps or crusts would also be prone to vulcanian explosions. In contrast, near-molten glowing surfaces with lava fountains and hornitos are favourable to strombolian eruptions, as were witnessed during 1984–2010.

Thus, Arenal appears to be rather a case of conduit fed at variable discharge rate in which persistent explosive activity could be controlled by changes in magma ascent conditions and thus also in decompression rates, as has been observed at volcanoes such as Tunguragua (Wright et al. 2012), Sakurajima (Miwa et al. 2009; 2013) or Galeras (Bain et al. 2021). Unfortunately, the low time resolution of the available discharge rates and magma ascent rates means that neither the short-term shifting/coexistence of vulcanian/strombolian explosivity, nor the shifting between effusive and explosive activity, can be constrained with any confidence. Further investigation is thus needed in order to better understand the evolution of the explosive activity at Arenal, which is a key aspect for volcanic hazard assessment.

Source models of explosion quakes and volcanic tremor

Many physical mechanisms have been proposed in the literature as possible sources of low-frequency seismicity at Arenal. The first model proposed to explain harmonic tremor was the acoustic resonance of a 1D conduit or organ pipe (Alvarado et al. 1988b). In this model, the frequencies of overtones are integer multiples of the fundamental mode frequency. However, this model suffers a weakness because it is difficult to produce the vibration of many overtones simultaneously. Likewise, the model of magma wagging (Bercovici et al. 2013), which describes the oscillations of a magma column surrounded by an annulus of bubbles, has a weakness because it can generate only one spectral peak and cannot produce the long-lasting oscillations observed at Arenal. Instead, Dmitrieva et al. (2013) developed a frictional faulting model to explain harmonic tremor with strong frequency gliding observed before eruptions of Redoubt volcano. In this model, the magma extrusion produces repeating stick–slip events in the shallow part of the conduit, with a recurrence interval decreasing with increasing stress rate (Denlinger and Moran 2014). However, this model cannot explain acoustic observations because its mechanism does not include a gas phase. Furthermore, it should produce a relationship between discharge rate and mean duration of tremor, but no correlation had been observed at Arenal.

Hagerty et al. (2000) and Julian (2000) proposed a mechanism to explain harmonic tremor at Arenal based on a nonlinear excitation by fluid flow through a conduit. This model is interesting because it can produce many regularly spaced spectral peaks without resonance, and it can generate the phenomena of period doubling and chaotic behaviour (Julian 1994). However, Rust et al. (2008) demonstrated that vibrations of the channel walls can occur if the velocity of the fluid is higher than a threshold, which depends on conduit dimensions and properties of the fluid and surrounding rocks. This condition can only be reached by high-velocity gas flow through a very long channel, which is unlikely to be the case at Arenal. Furthermore, the phenomenon of period doubling has never been observed in the hundreds of hours of tremor analysed by the authors of this paper, suggesting that this mechanism does not apply to Arenal.

Instead, we present a source model (Fig. 12) that explains most of the features of explosion quakes and tremor observed by seismic, acoustic and Doppler radar measurements, since the late 1990s (Lesage et al. 2006; Valade et al. 2012; Almendros et al. 2014). In this model, the viscous, degassed and cold cap, which closes the conduit, traps the exsolved gas and produces overpressure that can be released through fractures in the plug. If the release is sudden and involves a large amount of gas, the flow-induced vibrations produce an explosion quake (Fig. 12a and b). The gas flow associated with an explosion contains varying amounts of ash resulting from magma fragmentation or from plug rubble. It can produce either a radar signal (for ash-laden gas output), or not (for ash-free gas output) (Valade et al. 2012). Alternatively, the gas can escape intermittently through fractures that open when the pressure threshold is overcome and close when the pressure decreases. This process produces repetitive pressure pulses in the atmosphere that are audible as “chugs” (Melson 1989; Barboza and Melson 1990; Benoit and McNutt 1997). It also generates standing interface waves (Ferrazzini and Aki 1987) in the conduit that stabilise the period of pressure oscillations by a feedback mechanism (Lesage et al. 2006) (Fig. 12c and d). By analogy, this process that generates harmonic tremor has been called the “clarinet model”, in which the fractures act as a reed and the conduit is the resonating tube (Lesage et al. 2006). The repetitive seismic pulses generated by the opening and closing of fractures are the origin of the evenly spaced spectral peaks produced by a Dirac comb effect (Lesage et al. 2006). This implies that the fundamental frequency is the inverse of the repetition period and the spectrum is equal to the product of a Dirac comb in frequency by the spectrum of individual pulses. The more regular the repetition period of pulses, the larger the number of spectral peaks. The fracture is not strongly damaged for mild explosions and can still act as a valve, producing harmonic tremor in the coda of the event. If the fracture is blocked by rubble, the feedback mechanism is no longer efficient, the pressure pulses become irregular, and they generate spasmodic tremor. Transitions from spasmodic to harmonic can thus be explained by the evolving state of the fractures. The fundamental frequency of the conduit depends upon its length and on the velocity of the interface waves, which is related to the acoustic velocity of the magma. Small variations of the gas content in the magma can produce large changes in the velocity and thus frequency gliding (Sturton and Neuberg 2003).

Fig. 12
figure 12

source model for explosion quakes and tremor. a The pressure increases below the fractured, rigid cap. b When the pressure overcomes the strength of the fractures, gas can be released in two different manners. If the amount of gas flowing through the open fracture is large, the pressure drop induces magma fragmentation and varying quantities of pyroclasts are expelled during the explosion. Seismic waves are emitted from a few hundreds of metres below the crater. If the gas is released by repetitive pulses through the fracture, it generates standing waves in the conduit the period of which stabilises the pressure oscillations below the plug. The seismic waves are radiated at a very shallow depth by the fracture opening and closure and produce harmonic tremor by the Dirac comb effect (adapted from Fig. 8b of Valade et al. 2012). c The upper part of a conduit contains magma and gas bubbles. The solid cap in the crater and the exsolution level form closed ends of this pipe where standing tube waves can oscillate. d The existence of two conduits below the crater explains many observed features (e.g. double sets of spectral peaks with independent gliding, simultaneous explosion and tremor), modified from Lesage et al. (2006). e View of Arenal volcano from the NW. The fumaroles indicate the existence of two vents in the active crater C. The old crater D stands on the background of the picture. Photograph courtesy by Luis Madrigal (February 2004) and modified from Lesage et al. (2006). f Example of tremor showing two sets of spectral peaks with independent gliding behaviours modified from Lesage et al. (2006). Vertical component waveform in arbitrary units above and spectrogram below

Sketch of the plumbing system and

Full size image

The episodes of tremor containing two systems of overtones with independent frequency glidings (Fig. 12f) and the apparent lack of relationship between explosions and tremor could be explained by a double source associated to two magmatic conduits (Fig. 12d and e) (Lesage et al. 2006). If the two conduits are not linked, they form independent resonators (Lesage et al. 2006). Alternatively, if an acoustic link between the conduits is established (Fig. 12d), they can act as one resonator, with greater length and lower fundamental frequency.

The results of moment tensor inversions indicate that the source of explosions is slightly deeper below the crater (~ 200 m) than those of tremor (~ 100 m) (Davi et al. 2010, 2012). For both types of event, the position and depth of the sources are estimated through a probabilistic approach, using a grid search and taking the topography into account. The uncertainties are about ± 50 m in each direction (Davi et al. 2012), making the difference between locations of explosions and tremor significant. The seismic energy associated with tremor is mainly radiated at the top of the conduit (Jousset et al. 2003) and propagates as surface waves in the very heterogeneous shallow structure with rough topography. Thus, small modifications of the source position, or mechanism, may produce large variations in the path, back-azimuth and waveform. Moreover, two components of the wavefield which follow different travel paths may produce rapidly varying interferences as described by Almendros et al. (2012). However, seismic radiations associated with explosions occur at larger depth than the tremor and generate body waves, which propagate deeper in a more homogeneous medium. Their propagation is less sensitive to initial conditions at the source, and the histogram of their back-azimuth is more peaked than in the case of tremor (Almendros et al. 2014).

Acoustic observations: new insights into explosive activity at Arenal

A model of strombolian explosions can explain the waveforms of stacked small pressure transients. A good fit is obtained for the first positive pulse and for the second larger negative peak, which is a characteristic feature of this type of explosion, such as at Stromboli (Ripepe and Marchetti 2002; Ripepe et al. 2002) or Erta Ale (Bouche et al. 2010), and was modelled by Vergniolle and Brandeis (1996) (Fig. 11b). Furthermore, the maximum radiated acoustic pressure and gas fluxes for the five recorded explosions are much smaller at Arenal than for a pure vulcanian explosion, associated with a fully closed vent, such as at Sakurajima (e.g. Johnson and Miller 2014; Kim et al. 2015; De Angelis et al. 2019). This is in agreement with visual observations of the strongest explosion during our 2005 experiment, whose ash cloud was less than 1 km high. The small height of these gas-and-ash explosions suggests an open vent, such as for a strombolian explosion, rather than a closed vent as for a vulcanian one (Barboza and Melson 1990).

However, the seismic source models proposed for the explosions and tremor require a moderate pressurisation below a thin and weak solid cap. The contrasting view provided by acoustic and seismic records and models suggests a need to assess in the future whether a yet unknown mechanism could reconcile both interpretations, or if two different vents can operate in close proximity and at the same time with such a different pattern. A different pattern of explosions was also described at Stromboli between the southwest and northeast vents (Ripepe 1996) and often between the north and the south vents at Yasur (Iezzi et al. 2019). In addition, the Santiaguito open-vent system is well known to feed ash plumes and lava flow at the same time (since 1968) (Gottschämmer et al. 2021). However, in these cases, the differences in patterns were not as marked as between a fully open and a mechanically closed vent. The destruction of the lava cap during a vulcanian explosion is likely to produce a detectable seismic signal, whereas a strombolian explosion may or may not. Furthermore, no seismic signal is associated with the main infrasonic pulse of the strombolian explosions at Erta Ale (Bouche et al. 2010) and of some mildly explosive events at Stromboli (Vergniolle and Brandeis 1996), including puffing, because the radiation of the sound wave occurs mainly in the atmosphere. Therefore, the lack of seismicity associated with the small pressure transients at Arenal cannot be a decisive argument to choose between a mechanism based on a closed vent, likely to produce seismic events and a fully open vent.

Additionally, the existence of multiple vents and “hornitos” observed during the last stage of Arenal eruption could reconcile the two different types of acoustic signals and the source mechanisms inferred from the seismic, visual and Doppler radar observations (Lesage et al. 2006; Valade et al. 2012). While the small pressure transients could be generated by small strombolian explosions at some open vents accompanied by the build-up of hornitos, they could also result from an unknown nonlinear process associated with the closing and opening of a fracture, and not be driven by large bubbles as for purely strombolian activity. In this case, a relatively thin solid and fractured magma layer would induce moderate pressurisation of the underlying system at other vents and would play an important role in the mechanisms of the explosions and the generation of harmonic tremor. The behaviour of such a plumbing system, with simultaneously open and closed vents, probably connected at depth, would require further investigations.

In any case, the existence of the small acoustic waveforms, asymmetrical in favour of the rarefaction peak, which are typical of strombolian explosions, poses the question of the mechanisms of gas release at volcanoes with intermediate magma composition. Although the newly analysed acoustic measurements carried out in February 2005 only relate to the eruptive activity observed during that period, they provide a new insight that requires further studies and perhaps re-assessment of what was known so far about the eruptive mechanisms at Arenal, as well as for similar volcanoes.

Velocity models, site effects and deformation

Studies regarding Arenal edifice have provided general knowledge of its influence on the seismic wave field. At a global level, it was possible to adapt seismic antenna methods in order to locate the source of the seismic-volcanic signals with greater precision by considering site and topographic effects. However, there are some important aspects to be solved, such as: 1) the derivation of a deeper and more complete velocity model than the current one; 2) a better understanding, from a geotechnical perspective, of the behaviour of the layers that constitute the basement of the volcanic structure, and their response to the load exerted by the volcano, as well as its possible consequences; 3) monitoring of the volcano's slopes in order to know the long-term effect of the load exerted by the lava field and its stability; and 4) the necessity to understand deep processes related to magmatic plumbing systems that can be reconciled with geochemical knowledge.

Conclusions

We here review Arenal last eruptive phase (July 1968-October 2010), a phase spanning nearly 42 years of open-vent eruptive activity. This review considers literature and work for over 50 years of research at Arenal into geophysical aspects of the eruption. We have supported this review as well with a retrospective analysis of previously un-used data. Together, this collation of findings allowed us to:

i) Better understand the long-term persistent open-vent eruptive activity at Arenal volcano, especially during the declining phase. The eruption developed in seven stages that we characterise in this review in terms of the coupled eruptive-effusive processes, discharge rates, morphological changes and geophysical studies. Petrological, geochemical and rheological studies contributed to the understanding of a long-lasting open-vent activity in intermediate-composition (basaltic-andesitic) volcanism, particularly basaltic–andesitic, and to the modelling of lava flow dynamics. Deformation studies revealed the complexity of the deformation field in an open-vent volcano combining long-term subsidence with loading processes resulting from the near-constant effusion of lava flows and explosive activity, as well as the influence of erosion, reworking and bulking processes. The combination of seismic, Doppler radar and acoustic observations led to a new model of the source of tremor signals and their relation to explosive activity, as well as to reveal its complex nature.

(ii) Highlight outstanding aspects that remain poorly understood at Arenal such as:

  1. a)

    the alternation and/or coexistence of strombolian and vulcanian explosive regimes coupled with an almost persistent effusive activity,

  2. b)

    the nature of the lava hosted inside the active crater,

  3. c)

    the magma ascent processes as well as the configuration of the magmatic feeding system and structure of the volcanic edifice.

(iii) Offer new insights into future research topics, some of which could be addressed by the detailed study of the products of the 1968–2010 eruption and others that could be addressed during the next eruption. Efforts towards a better description of the explosive products, together with seismic and acoustic observations and modelling, are important for understanding the dynamics in the Arenal volcano conduit and thus for understanding other basaltic-andesitic systems. Geophysical studies aimed at modelling the structure of the volcanic edifice and the feeder system would be fundamental for a better understanding of unrest processes and monitoring before and during a future eruptive phase.

Arenal's 1968–2010 open-vent eruption contributed greatly to the development of modern volcano-monitoring and research in Costa Rica. During the eruption, substantial improvement in different monitoring methods at open-vent volcanoes has occurred in Costa Rica and elsewhere, focusing on the construction of continuous and comprehensive databases that allow adequate tracking of activity and interpretation of processes driving that activity as well as any changes therein. Forecasting of eruptions at the Costa Rican volcanoes of Turrialba, Poás, and Rincón de la Vieja volcanoes during the last two decades has been proof of how the 1986–2010 eruption of Arenal greatly advanced the scientific level of Costa Rican observatories so as to arrive at the current status. These advances, and solutions to open questions, will now allow observatories to be better prepared for future eruptive activity at Arenal and similar open-vent systems elsewhere.