A Modiﬁed Sandpile Model for Simulating Lava Fountains at Mt Etna

Giuseppe Nunnari

Dipartimento di Ingegneria Elettrica, Elettronica e Informatica, Università degli Studi di Catania,

Viale A. Doria, 6, 95125 Catania, Italy

Keywords:

Sandpile Model, Self-Organized Criticality (SOC), Power-Law Distribution, Volcanic Lava Fountains,

Inter-Event Time Distribution.

Abstract:

This study aims to achieve two primary objectives. Firstly, it offers empirical evidence on the statistical dis-

tribution of inter-event times for lava fountains at Mt. Etna between 2011 and 2022, revealing that these times

follow a power-law distribution, which supports the hypothesis that volcanic energy release exhibits dynam-

ics characteristic of Self-Organized Criticality (SOC) systems. Secondly, it introduces a modiﬁed version of

the classic Bak-Tang-Wiesenfeld (BTW) model, speciﬁcally adapted to simulate the inter-event times of lava

fountains in volcanic environments like Mt. Etna. Although the proposed model is straightforward and offers

initial insights, it remains preliminary. Further development is needed to enhance its accuracy and extend its

applicability to more complex volcanic systems.

1 INTRODUCTION

Many natural systems exhibit self-organization at the

edge of phase transitions, a phenomenon known as

self-tuned phase transitions. Such systems have been

observed across various domains, including labora-

tory fusion plasmas, solar physics, and magneto-

spheric physics (Watkins et al., 2015). In geophysics,

many systems display SOC characteristics (Corral

and Gonzales, 2019). It is now widely recognized

in seismology that seismic energy release follows a

power-law distribution, with different regions exhibit-

ing distinct representations of the Gutenberg-Richter

law of earthquake magnitude (Corral and Gonzales,

2019). The presence of power laws is a strong in-

dicator of Self-Organized Criticality (SOC) (Jensen,

1998), suggesting that seismic energy release is a

classic example of a SOC system.

Given the SOC characteristics observed in seis-

mic energy release, it is plausible to hypothesize that

volcanic energy release might follow similar dynam-

ics. Volcanic activity encompasses a wide range of

eruption styles, with signiﬁcant variations in the vol-

ume of ejected material and the intervals between

eruptions. Previous studies have identiﬁed power-

law relationships between the Volcanic Explosivity

Index (VEI) and eruption frequency, drawing a par-

allel to the Gutenberg-Richter law observed in earth-

https://orcid.org/0000-0002-7117-3174

quake magnitude distributions (Butters et al., 2017).

This suggests a potential SOC behavior in volcanic

activity as well.

This study speciﬁcally investigates the inter-event

times of lava fountains at Mt. Etna, aiming to demon-

strate their adherence to a power-law distribution.

Mt. Etna is one of the most active volcanoes in the

world, making it an ideal candidate for studying vol-

canic energy release through lava fountains. To model

this phenomenon, we adapt the Sandpile model (Bak

et al., 1987), a cellular automaton originally devel-

oped to illustrate SOC, which has been widely used to

model critical phenomena in complex systems. Cel-

lular automata, with their discretized time, space, and

internal states, can produce complex structures de-

spite their simple local evolution rules. They have

been applied to model non-equilibrium phenomena

such as crystal growth, ﬂuid dynamics, and transport

processes governed by the Navier-Stokes equations.

In the context of volcanic systems, cellular automata

have been employed to simulate intricate phenomena,

including lava ﬂows (Vicari et al., 2007) and magma

ascent (Piegari et al., 2011). In this study, we explore

an adapted version of the BTW model to interpret the

inter-event times of Lava Fountains (LF) at Etna, pro-

viding further experimental evidence that volcanic en-

ergy release may follow power-law distributions.

700

Nunnari, G.

A Modiﬁed Sandpile Model for Simulating Lava Fountains at Mt Etna.

DOI: 10.5220/0013043500003822

Paper published under CC license (CC BY-NC-ND 4.0)

In Proceedings of the 21st International Conference on Informatics in Control, Automation and Robotics (ICINCO 2024) - Volume 1, pages 700-705

ISBN: 978-989-758-717-7; ISSN: 2184-2809

2 LAVA FOUNTAIN

PHENOMENON

Lava fountains, such as those observed at Mt Etna,

can be understood as resulting from the gradual ac-

cumulation of stress within the volcanic system. This

stress eventually triggers the ascent of molten mate-

rial, with gases playing a crucial role in facilitating

this rise. However, the mechanisms that control such

events are not fully understood (La Spina et al., 2017),

and their modeling using partial differential equations

(PDEs) presents signiﬁcant challenges.

Firstly, our understanding of the physical proper-

ties within a volcanic system is limited. The subsur-

face conditions, such as the structure and composition

of the magma chamber and conduit system, are poorly

constrained by available observational data, making it

difﬁcult to accurately parameterize these features in

mathematical models.

Secondly, magma itself is not a simple ﬂuid. It

exhibits non-Newtonian behavior, where its viscos-

ity varies depending on factors like temperature and

crystal content. The presence of solid particles fur-

ther complicates its ﬂow dynamics, introducing chal-

lenges that standard ﬂuid dynamics equations struggle

to capture accurately.

Furthermore, the interaction between rising

magma and the surrounding rock, coupled with the

complex degassing processes, introduces additional

layers of complexity. These interactions can affect

the magma’s ascent rate, pressure buildup, and ulti-

mately, the explosivity of eruptions. The exsolution

of gases from magma plays a crucial role in driving

ascent and eruption processes, but these processes are

highly non-linear and not fully understood.

Additionally, the stochastic nature of volcanic ac-

tivity, inﬂuenced by numerous random factors, com-

plicates the prediction of precise behaviors using de-

terministic models. This inherent variability, cou-

pled with the multitude of interacting processes, often

renders traditional PDE-based models inadequate for

capturing the full range of volcanic behaviors.

Given these challenges, alternative modeling ap-

proaches like the Sandpile model (Bak et al., 1987),

which is adapted in this study, offer valuable insights

into the statistical properties of volcanic phenomena.

This model allows for the exploration of eruption dy-

namics without requiring detailed physical parame-

ters, making it a powerful tool for studying the dis-

tribution of time intervals between eruptions.

3 DATA AND METHOD

3.1 Experimental Data

The dataset for this study comprises the start and

end dates of 118 lava fountain episodes at Mt Etna,

recorded between 2011 and 2022 (Calvari and Nun-

nari, 2022). These episodes have been meticulously

documented, offering a detailed chronicle of volcanic

activity during this period.

To analyze the temporal distribution of these lava

fountain episodes, we calculated the inter-event times,

deﬁned as the difference in days between the start of

one lava fountain and the start of the next. This metric

quantiﬁes the intervals between successive volcanic

activities, providing insights into the underlying tem-

poral patterns (Calvari and Nunnari, 2022).

Figure 1: Histogram of LFs inter-event times. The x-axis

represents the inter-event time in days, while the y-axis rep-

resents the frequency of occurrences. This histogram illus-

trates the distribution of intervals between consecutive lava

fountain episodes at Mt Etna from 2011 to 2022, highlight-

ing the variability and periodicity of volcanic activity.

Figure 1 presents the histogram of LFs inter-event

times, offering a visual representation of the temporal

distribution of lava fountain episodes. This histogram

reveals the variability and periodicity in the volcanic

activity at Mt Etna over the studied period, which is

crucial for understanding the system’s dynamics.

3.2 Some General Features of SOC

Systems

Self-organized criticality (SOC) is a property of dy-

namical systems that naturally evolve to a critical

state, where a minor event can lead to signiﬁcant

consequences. This concept has been applied to a

wide range of physical, biological, and social sys-

A Modiﬁed Sandpile Model for Simulating Lava Fountains at Mt Etna

701

tems. Here are some general features of SOC sys-

tems:

• Scale Invariance. SOC systems exhibit behav-

ior that is independent of the scale at which they

are observed. This means that similar patterns or

structures can be found at different magnitudes,

indicating a fractal-like nature. For a function of

the type f (x) = x

, the relative variation

f (kx)

f (x)

does not depend on x.

• Power-law Distributions. The size and fre-

quency of events in SOC systems typically fol-

low a power-law distribution. This implies that

small events are exceedingly common, while large

events are rare but possible. One method to ex-

perimentally determine if a phenomenon could be

generated by an SOC system is to measure a char-

acteristic of the phenomenon, in space or time,

and check if its probability distribution follows a

power-law. However, even if a power-law distri-

bution is observed, this does not constitute rigor-

ous proof of SOC behavior, as the absence of a

characteristic scale might be obscured or not evi-

dent in the dataset. This means that while observ-

ing a power-law distribution may suggest SOC be-

havior, it is not deﬁnitive evidence of it; other un-

derlying factors or data limitations could obscure

the true nature of the system.

• Threshold Dynamics. SOC systems often have a

critical threshold that, when reached, can trigger

a chain reaction or cascading events. The system

naturally evolves towards this critical point with-

out external tuning.

• Long-range Correlations. In SOC systems, dis-

tant parts of the system can become correlated,

meaning that a change in one part of the system

can inﬂuence another part, regardless of the phys-

ical distance.

• Intermittent Activity. The activity in SOC sys-

tems is typically sporadic, with periods of relative

calm interspersed with bursts of intense activity.

This intermittent behavior is a hallmark of critical

systems.

Understanding these features helps in identifying

and analyzing SOC behavior in various complex sys-

tems, from earthquakes to ﬁnancial markets.

3.2.1 The BTW Model for Simulating SOC

Systems

The BTW (Bak-Tang-Wiesenfeld) Sandpile model is

a mathematical framework commonly used to sim-

ulate self-organized criticality (SOC) across various

systems, including volcanic activity. The model oper-

ates on a grid, where each cell can hold a certain num-

ber of grains. When a cell exceeds a critical threshold,

it "topples," distributing grains to neighboring cells

and potentially causing a chain reaction of toppling

events, known as an avalanche.

In its standard form, the Sandpile model operates

as follows:

1. Initialization. A square grid of size pile_width ×

pile_width is initialized. Each cell in the grid is

randomly assigned up to three grains.

2. Grain Addition. Grains are added to the grid one

at a time at random positions. If adding a grain

causes a cell to exceed three grains, an avalanche

is triggered.

3. Avalanche Process. During an avalanche, four

grains are removed from the overloaded cell: one

grain is added to each of its four neighboring cells

(top, right, bottom, left). If any of these neigh-

boring cells exceed three grains as a result, they

also topple, continuing the avalanche until no cell

exceeds three grains.

4. Boundary Conditions. Grains that topple off the

edge of the grid are lost and not considered in sub-

sequent steps.

5. Simulation Loop. The process of adding grains

and resolving avalanches continues until the spec-

iﬁed number of grains has been added to the grid.

3.2.2 Adaptation of the BTW Model for

Simulating Lava Fountains

To better simulate lava fountains, the adapted BTW

model introduces critical modiﬁcations to the bound-

ary conditions and the resistance matrix, which are

necessary to capture the unique behaviors observed

in volcanic activity. Below, we outline the major dif-

ferences between the standard BTW model and our

adapted version.

Consider a 2D grid of size N ×N. Each cell in the

grid is denoted by (i, j) with i, j ∈ {1, 2, . . . , N}. The

rules governing the accumulation and release of stress

are as follows:

1. Stress Accumulation Rule. Stress, representing

the buildup of volcanic pressure, is incrementally

added to a random position in the bottom row dur-

ing each simulation step:

S(N, j) → S(N, j) + 1,

where j is randomly chosen from {1, 2, . . . , N}.

2. Resistance Matrix.

– The resistance matrix R(i, j) decreases from the

bottom to the top of the grid.

ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics

702

– The resistance on the lateral boundaries (left and

right edges) is set to inﬁnity:

R(i, j ) =

(

∞, if j = 1 or j = N

min

i−1

N−1

× (R

max

− R

min

), otherwise

A visual representation of the resistance matrix is

shown in Figure 2.

Figure 2: An example of the resistance matrix adapted for

simulating volcanic activity, illustrating how resistance de-

creases with height and inﬂuences the dynamics of stress

accumulation and release.

3. Stress Redistribution Rule. When the stress

S(i, j) in a cell exceeds the resistance R(i, j), an

event occurs:

S(i, j) → S(i, j) − 4,

and the excess stress is redistributed to the four

neighboring cells:

S(i + 1, j) → S(i + 1, j) + 1,

S(i − 1, j) → S(i − 1, j) + 1,

S(i, j + 1) → S(i, j +1) + 1,

S(i, j − 1) → S(i, j −1) + 1.

This redistribution of stress can trigger further

events in neighbouring cells, leading to a chain re-

action, referred to as an avalanche.

4. Boundary Conditions. The stress in cells at the

top and bottom rows is lost if it attempts to redis-

tribute outside the grid:

S(i, j) → 0 for i = 1 or i = N.

5. Stop Condition. The simulation is terminated

when the accumulated stress in the entire grid

reaches a threshold value.

These modiﬁcations, particularly the introduction

of a resistance matrix that decreases with height and

the speciﬁc boundary conditions, are critical for repli-

cating the self-organized criticality observed in vol-

canic systems. The next subsection will explore the

implications of these changes for understanding the

dynamics of lava fountains at Mt Etna.

4 RESULTS

This section presents the analysis of lava fountain

(LF) inter-event times and their compatibility with

a power-law model. Section 4.1 evaluates whether

the observed inter-event times of LFs at Mt. Etna ﬁt

a power-law distribution, providing the correspond-

ing model parameters. Section 4.2 examines whether

the inter-event times from simulated LFs, using the

adapted sandpile model, exhibit similar power-law

behavior, thereby assessing the model’s ability to

replicate the characteristics of the actual volcanic

events.

4.1 A Power-Law for the Mt Etna LF

Inter-Event Times

To determine if the lava fountains at Mt. Etna exhibit

features characteristic of a Self-Organized Criticality

(SOC) system, we analyzed the inter-event times of

lava fountains recorded from 2011 to 2022. Following

the analysis by (Calvari and Nunnari, 2022), the inter-

event times were ﬁtted to a power-law distribution,

resulting in the parameters α = 1.72, x

min

= 5.39,

= 0.19, and σ

min

= 6.49. The p-value of 0.46 in-

dicates that the power-law model is statistically plau-

sible, suggesting that the inter-event times lack a char-

acteristic scale, which is consistent with the behavior

of SOC systems.

Figure 3 shows the empirical complementary cu-

mulative distribution function (CCDF) for the LF

inter-event times and the corresponding power-law ﬁt,

highlighting the ﬁt’s alignment with the data.

4.2 Simulating the Adapted Sandpile

Model

To assess whether the adapted sandpile model can

replicate the inter-event times of lava fountains ob-

served at Mt. Etna, we conducted simulations over

50,000 cycles on a 64x64 grid. We tracked both the

number of avalanches and the grains expelled from

the top central portion of the grid, which we identiﬁed

as simulated lava fountain events. The times between

these events are considered as the inter-event times of

simulated lava fountains.

A Modiﬁed Sandpile Model for Simulating Lava Fountains at Mt Etna

703

Figure 3: Empirical complementary cumulative distribution

function (CCDF) (blue circles) and its power-law model ﬁt

(dotted black line) for the inter-event times of lava fountains

(LF) at Mt. Etna between January 2011 and July 2022. The

symbol ’*’ indicates the x

min

for the considered dataset.

4.2.1 Simulating the Avalanche Size

Figure 4 illustrates the distribution of observed

avalanches by size as recorded by the model. This

distribution is a key indicator of the model’s behavior

and its potential alignment with SOC characteristics.

Figure 4: Number of observed avalanches versus their size

for the adapted 64x64 grid size sandpile model.

The empirical CCDF for avalanche sizes, along

with the ﬁtted power-law model, is presented in Fig-

ure 5. The parameters obtained (α = 1.70, x

min

= 8)

and the p-value of 0.60 suggest that the model re-

produces the avalanche dynamics characteristic of

the original Bak-Tang-Wiesenfeld (BTW) sandpile

model, reinforcing the model’s validity.

4.2.2 Simulating the LF Inter-Event Times

Over 50,000 simulation cycles, the model recorded

68 lava fountain events, allowing for the analysis of

Figure 5: Empirical CCDF of avalanche sizes and the cor-

responding power-law ﬁtting model.

inter-event times. The power-law ﬁt for the simulated

LF inter-event times yielded parameters of α = 1.59,

min

= 47, σ

= 0.25, and σ

min

= 180.64. The p-

value of 0.09 suggests marginal plausibility, indicat-

ing that while the model approximates the observed

distribution, it may not fully replicate the statistical

characteristics of the actual LF inter-event times at

Mt. Etna. This result suggests limitations in the

model’s ability to fully capture the dynamics of the

volcanic system.

Figure 6: CCDF of the simulated LF inter-event times using

the adapted 64x64 grid size sandpile model.

5 CONCLUSIONS

This study demonstrates that the proposed sandpile

model effectively reproduces the avalanche dynam-

ics characteristic of the Bak-Tang-Wiesenfeld (BTW)

model, yet it only partially succeeds in simulating the

inter-event times of the lava fountain phenomena ob-

ICINCO 2024 - 21st International Conference on Informatics in Control, Automation and Robotics

704

served at Mt. Etna. Several factors contribute to the

model’s partial success in replicating real-world lava

fountain dynamics:

Firstly, the simpliﬁed stress propagation mecha-

nism in the model likely fails to fully capture the com-

plex dynamics inherent in real lava fountains. This

simpliﬁcation may overlook critical interactions and

feedback processes. Future reﬁnements could include

integrating more sophisticated stress redistribution al-

gorithms or exploring alternative stress propagation

models that better reﬂect the complexity of volcanic

processes.

Secondly, the limited grid size (64x64) and the

number of simulation cycles (50,000) may have con-

strained the model’s ability to accurately simulate the

complexity of lava fountains. Larger grids and ex-

tended simulation periods could capture a broader

range of inter-event times, potentially leading to a

better approximation of real-world dynamics. Future

studies should consider scaling up these parameters to

assess their impact on the accuracy of the simulation

results.

Furthermore, the comparison with the more com-

plex model by (Piegari et al., 2011) underscores the

inherent trade-offs between model simplicity and ac-

curacy. While our model offers a tractable frame-

work for initial investigations, it may lack the detailed

mechanisms necessary for higher ﬁdelity simulations.

This balance between simplicity and detail highlights

the need for ongoing model reﬁnement and validation

to better capture the nuances of volcanic activity.

Overall, this work represents a valuable initial

attempt to simulate a complex geophysical phe-

nomenon using a simpliﬁed model. The insights

gained from this study offer a foundation for future

research, which could focus on reﬁning the model to

better capture the nuances of lava fountain dynamics

and exploring alternative approaches to improve sim-

ulation accuracy. Continued development and valida-

tion of the model will be crucial in enhancing its abil-

ity to reproduce real-world lava fountain phenomena

and contribute to our understanding of these complex

systems.

Future research directions may include:

• Further reﬁning the stress propagation mechanism

to better mirror the dynamics observed in actual

lava fountains, possibly by incorporating more

complex algorithms or alternative modeling ap-

proaches.

• Systematically testing the model with larger grid

sizes and extended simulation periods to explore

how these factors inﬂuence the ﬁdelity of the sim-

ulated lava fountain dynamics, potentially leading

to more accurate and reliable predictions.

• Investigating additional external factors, includ-

ing environmental conditions and geological vari-

ations, that could signiﬁcantly inﬂuence lava

fountain behavior, offering a more holistic under-

standing of the processes involved.

ACKNOWLEDGEMENTS

We would like to thank the INGV-OE scientists and

technicians for the monitoring network maintenance,

and especially Dr Sonia Calvari for providing essen-

tial information for this work.

REFERENCES

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-

organized criticality: an explanation of 1/f noise.

Physical Review Letters, 59(4):381–384.

Butters, O., Sarson, G., and Bushby, P. (2017). Effects

of magma-induced stress within a cellular automa-

ton model of volcanism. Journal of Volcanology and

Geothermal Research, 341:94–103.

Calvari, S. and Nunnari, G. (2022). Etna output rate dur-

ing the last decade (2011–2022): Insights for hazard

assessment. Remote Sensing, 14:1–16.

Corral, A. and Gonzales, A. (2019). Power law size dis-

tributions in geoscience revisited. Earth and Space

Science, 6:673–697.

Jensen, H. (1998). Self-Organized Criticality: Emergent

Complex Behavior in Physical and Biological Sys-

tems. Cambridge University Press, Cambridge, UK.

La Spina, G., dè Michieli Vitturi, M., and Clarke, A.

(2017). Transient numerical model of magma ascent

dynamics: application to the explosive eruptions at the

soufrière hills volcano. Journal of Volcanology and

Geothermal Research, 336:118–139.

Piegari, E., Di Maio, R., Scandone, R., and Milano, L.

(2011). A cellular automaton model for magma as-

cent: Degassing and styles of volcanic eruptions.

Journal of Volcanological and Geothermal Research,

30:22–28.

Vicari, A., Herault, A., Del Negro, C., Coltelli, M.,

Marsella, M., and Proietti, C. (2007). Modelling of

the 2001 lava ﬂow at etna volcano by a cellular au-

tomata approach. Environmental Modelling and Soft-

ware, 22:1465–1471.

Watkins, N., Pruessner, G., S.C., C., Crosby, N., and H.J., J.

(2015). 25 years of self-organized criticality: concepts

and controversies. Science Review.

A Modiﬁed Sandpile Model for Simulating Lava Fountains at Mt Etna

705