# Transition from volcano-sagging to volcano-spreading

Eoghan P. Holohan, Sam Poppe, Audray Delcamp, Paul K. Byrne, Thomas R.

Walter, Benjamin van Wyk de Vries, Matthieu Kervyn

► To cite this version:

Eoghan P. Holohan, Sam Poppe, Audray Delcamp, Paul K. Byrne, Thomas R. Walter, et al. Transition from volcano-sagging to volcano-spreading. Earth and Planetary Science Letters, 2023, 604, 10.1016/j.epsl.2023.118012. insu-04197048

HAL Id: insu-04197048

https://insu.hal.science/insu-04197048

Submitted on 6 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Transition from volcano-sagging to volcano-spreading

Eoghan P. Holohan<sup>a,*</sup>, Sam Poppe<sup>b</sup>, Audray Delcamp<sup>c</sup>, Paul K. Byrne<sup>d</sup>, Thomas R. Walter<sup>e</sup>,
Benjamin van Wyk de Vries<sup>f</sup>, Matthieu Kervyn<sup>c</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Article history:**

Received 26 August 2022

Received in revised form 2 January 2023

Accepted 11 January 2023

Available online 1 February 2023

Editor: C.M. Petrone

**Keywords:**

volcano deformation

volcanic rift zones

analogue models

Digital Image Correlation

Mars

Earth

## ABSTRACT

Many volcanoes exhibit evidence of long-term deformation driven by gravity. This deformation influences the development of magmatic plumbing systems and volcanic vent locations, and it manifests as earthquakes and ground surface motions. Two end-member deformation styles are volcano sagging (flexure) and volcano spreading. Previous modeling studies indicated that the mechanical and geometric properties of basement rocks below the volcano can control whether a volcano spreads or sags, but these works differed on exactly how. Furthermore, evidence in nature that a volcano may progress from sagging to spreading as it grows has not been tested experimentally. Using Digital Image Correlation (DIC) and a revised scaling approach, we here revisit the classic experiments in which a sand-plaster cone is emplaced upon a basement comprising an upper brittle layer and a lower ductile layer. In a first experiment series, we emplaced the cone instantaneously. By normalizing both brittle and ductile layer thicknesses to cone height, we provide new insight into the dimensionless geometric controls on spreading and sagging. In a novel second experiment series, we emplaced the cone incrementally. We show that as a growing cone migrates through the dimensionless geometric parameter space, its deformation style changes successively in time from sagging to spreading. DIC shows that the horizontal velocity of the cone flank increases linearly or non-linearly during cone construction, but it decreases exponentially after construction ceases. Despite their geometric simplification and the omission of several sedimentological, magmatic, and regional-tectonic processes, our models' results are compatible with a range of observations of volcanoes on Earth and Mars. In particular, they help to explain: (i) why lithosphere-scale loading results in sagging rather than spreading, and (ii) why ocean island volcanoes of >2-3 km in height develop stellate forms and rift zones.

© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The deformation of volcanoes is closely intertwined with their growth via magma intrusion and eruption. Forces driving volcano deformation stem from a combination of magmatic activity, hydrothermal activity, regional tectonics, and gravitational instability (Chiocci et al., 2011; Moore, 1970; Neri et al., 2009; Poland et al., 2017). Such deformation occurs on both long (millennial to million-year) and short (daily to yearly) timescales. Long-term deformation of volcanic edifices is expressed primarily in volcano morphology and geological structure (Fig. 1). Short-term deformation

is now routinely documented on Earth by geodetic and seismic observations. Understanding the controls on, and the structural and dynamic consequences of, long-term volcano deformation is more challenging, but provides a basis for interpretation of short-term deformation and hence for hazard assessment of active volcanic systems worldwide (Poland et al., 2017).

Gravitational instability gives rise to two end-member styles of long-term deformation that are the subjects of our study: volcano sagging and volcano spreading (Borgia et al., 2000; McGovern et al., 2015). These styles are characterized by deformation of both the volcanic edifice and its underlying basement, and they are driven by gravitational loading from the edifice. Both deformation styles are characterized by low strain rates in the long term (10⁴–10⁶ years), although in the short term they can manifest

<footnote><sup>a</sup> UCD School of Earth Sciences, University College Dublin, Dublin 4, Ireland</footnote>

<footnote><sup>b</sup> Centrum Badań Kosmicznych Polskiej Akademii Nauk (CBK PAN), Bartycka 18A, 00-716 Warsaw, Poland</footnote>

<footnote><sup>c</sup> Department of Geography, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium</footnote>

<footnote><sup>d</sup> Department of Earth and Planetary Sciences, Washington University, Campus Box 1169, 1 Brookings Dr, Saint Louis, MO 63130-4899, United States of America</footnote>

<footnote><sup>e</sup> GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany</footnote>

<footnote><sup>f</sup> Laboratoire Magmas et Volcans, 5 rue Kessler, 63038 Clermont-Ferrand, France</footnote>

<footnote>* Corresponding author.
E-mail address: eoghan.holohan@ucd.ie (E.P. Holohan).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Volcanoes on Earth and Mars affected by long-term, gravity-driven deformation. Digital elevation models for Earth are from GMRT; those from Mars are from MOLA. These edifices display morphological or structural features interpreted to have been generated by spreading or sagging (flexure), or by some combination of spreading and sagging (Borgia and van Wyk de Vries, 2003; Byrne et al., 2013; Kervyn et al., 2014; McGovern et al., 2015; Platz et al., 2011).</figcaption></figure>

as higher-strain-rate events such as volcanotectonic or slow-slip earthquakes (Brooks et al., 2008; McGovern, 2007). Long-term deformation may also provide the framework for shorter-term mass movements, such as rotational landslides or debris avalanches (Moore et al., 1989), and it may either facilitate or inhibit the intrusion and eruption of magma (McGovern et al., 2015; Poland et al., 2017).

For the purposes of this study, the volcano “basement” may constitute anything from upper crustal materials immediately below an edifice to the entire lithosphere–asthenosphere mechanical system. On smaller scales, volcano sagging or spreading is commonly linked with weak upper-crustal strata, such as those containing clay, salt, or hydrothermally altered materials (Borgia et al., 2000). Such strata are located either at the volcano base or at some level below that, and they provide a detachment zone (or ‘décollement’) for relative motion between the volcano and stronger, more rigid crustal materials farther below (Byrne et al., 2013; Delaney et al., 1998). On the largest scale, volcano sagging results from the flexural-isostatic response of planetary lithospheres to loading by large volcanic edifices, and it is accommodated by ductile flow of the underlying asthenosphere (McGovern and Solomon, 1993; Walcott, 1970; Zucca et al., 1982).

Much of our understanding of gravitationally driven volcano deformation stems from analogue or numerical models (Borgia, 1994; Byrne et al., 2013; McGovern and Solomon, 1993; Merle and Borgia, 1996; Norini and Acocella, 2011; van Wyk de Vries and Matela, 1998). End-member volcano spreading typically produces near radially orientated faults with normal to oblique-normal or strike-slip motions that accommodate extension of the volcano flanks (Fig. 1). Ideally, with symmetrical spreading, these faults define a ‘flower-shaped’ pattern of near-radial, petal-shaped grabens and wedge-shaped horsts. Radial contraction commonly occurs around

the cone base, where it produces concentrically orientated folds and reverse (thrust) faults, as well as near-radially orientated conjugate strike-slip faults.

End-member volcano sagging generates the opposite spatial arrangement of extension and contraction. The volcano flanks undergo contraction and develop near-concentrically orientated faults with reverse motions. Ideally, these faults define flank terraces arranged in an imbricated or ‘fish-scale’ pattern (Fig. 1). Folds or reverse faults may also occur in a near-radial orientation to accommodate contraction, especially on the lowermost flanks. At and surrounding the volcano base is a flexural trough, outward from which lies a concentric zone of normal faulting and a flexural bulge.

The morphological and structural characteristics of these deformation styles occur at many volcanoes (Fig. 1). Long-term volcano spreading has been noted in the flank motions, rift zones, and basal structures of Hawaii, Etna, and Mt Cameroon volcanoes (Kervyn et al., 2014; Poland et al., 2017), although buttressing effects and edifice elongation cause deviation of the structures and motions from the ideal case. Spreading is suspected to cause the stellate planform morphology of many large volcanic seamounts and continental arc volcanoes on Earth (Delcamp et al., 2008; Mitchell, 2001), as well as the distinctive morphology of Tharsis Tholus on Mars (Platz et al., 2011). Long-term volcano sagging has been documented in the down-warped sea-floor and crust surrounding oceanic seamounts, such as at Hawaii (Zucca et al., 1982) and Tenerife (Watts et al., 1997). Other evidence of volcano sagging includes radial folds or wrinkle ridges, such as at the base of Concepción volcano, and flank terraces, such as seen at the Galapagos Platform on Earth and on Ascreaus Mons on Mars. Volcanoes such as Concepcion, Hawaii, and Olympus Mons display evidence of both spreading and sagging structures (Byrne et al., 2013), or,
<figure><img src="image_5.png" /><figcaption>(B) 'Incremental' Cone Growth</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Sketches of experimental set-ups for (A) instantaneous cone growth or (B) incremental cone growth. Note that (B) excludes syn-growth deformation, and only represents the volume addition during the cone growth experiments. The inset shows cone volume with time for a representative experiment (S-Sinc04).</figcaption></figure>

as in the case of Vesuvius, a possible hybrid geodetic behavior (De Matteo et al., 2022).

Analysis of structures within the annular fold belt that surrounds Concepción volcano suggests that the volcano initially underwent sagging followed by spreading (Borgia and van Wyk de Vries, 2003). Moreover, small oceanic seamounts typically have circular planform shapes, whereas large seamounts above a critical height of 2–3 km commonly display stellate planform shapes associated with volcanic rifts (Mitchell, 2001). These observations have led to the hypothesis that volcanoes could potentially transition from one deformation style to another with time as they grow (Borgia, 1994; Vogt and Smoot, 1984).

Past experimental and numerical modeling studies have proposed varying geometric thresholds at which spreading transitions to sagging (Borgia et al., 2000; Merle and Borgia, 1996; van Wyk de Vries and Matela, 1998). In this paper, we provide improved experimental constraints based on a revised scaling approach and a far greater number of experiments than in previous works to gain new insights into how geometric factors control whether a volcano spreads or sags. We also use Digital Image Correlation to quantify our experimental outcomes and to characterize surface displacements, strains, and velocities. Finally, we report new experiments simulating incremental volcano construction. These provide the first explicit demonstration of how a transition from basement-controlled sagging to spreading can occur during volcano growth.

## 2. Experimental methodology

### 2.1. Model set-up

Following the approach of Merle and Borgia (1996), we consider the natural system as a frictional-plastic cone of radius R and height H emplaced upon a basement comprising an upper frictional-plastic layer of thickness B and a lower Newtonian viscous layer of thickness D (Fig. 2 and Table 1). For simplicity, we neglected complexities arising from major magmatic intrusions and related hydrothermal systems, although we recognize that in nature these may have interactions and feedbacks with basement-

induced deformations (Le Corvec and Walter, 2009; Norini and Acocella, 2011). We also neglected the superposition of deformational events rooted at multiple levels in the basement (Byrne et al., 2013), effects of non-axisymmetric edifice shape (Kervyn et al., 2014), effects of lateral cone confinement ('buttressing') (Merle and Borgia, 1996) and the interaction with pre-existing or coevally developing regional tectonic structures (Grosse et al., 2020). The Discussion section briefly addresses some of these limitations.

The frictional-plastic ('brittle') material was a sand and plaster mixture at 90:10 by weight; this follows a Mohr-Coulomb failure criterion with a cohesion of ~60 Pa and a friction coefficient of ~0.8 (Poppe et al., 2021). The Newtonian viscous ('ductile') material was a polydimethylsiloxane (PDMS) (Weijermars, 1986). The PDMS was recycled due to the high number of model runs, such that it absorbed some plaster and fine sand particles. Dynamic viscosity measurements for 'clean' and 'dirty' PDMS remained within the same order of magnitude (4–8 × 10⁴ Pa), however. These materials were emplaced within a 50 cm-diameter ring as horizontal layers (Fig. 2), upon which a cone of height 5 cm and radius 8.7 cm was constructed.

We conducted two sets of experiments. In a first set of 75 experiments, the sand-and-plaster cone (of volume ~400 cm³) was constructed instantaneously at the angle of repose (30 ± 3°) by pouring (Fig. 2A). Experiment durations depended on deformation rate and varied from 120 to 300 minutes. In a second set of nine experiments, the cone was constructed incrementally at the angle of repose. Eight pours of volume ~50 cm³ at 45-minute intervals were made typically upon an initial seed cone of volume 3 cm³ (Fig. 2B).

For both experiment sets, the thicknesses of the ductile and brittle substrata were systematically modified, and some configurations were repeated to ensure reproducibility of results. To help visualize the surficial strains in the experiments, a plaster veneer (<1 mm thick) was sieved on top of each model. Additionally, several incremental growth models were wetted and quickly cross-sectioned at the end of certain stages of growth to monitor changes in the average thickness of the substrata during growth.
<caption>**Table 1**</caption>

Governing parameters and dimensionless Π-numbers for basement-controlled, gravity-driven deformation of a volcano. The parameter ranges of H, R, B, and D for natural volcanoes are taken from references cited in supplementary Table S6. Sources for density, cohesion, friction and viscosity estimates in nature are given also in the supplementary material. For the Reynolds number, minimum and maximum horizontal velocities used are 6 × 10⁻⁸ – 1 × 10⁻⁵ m s⁻¹ in experiment (from DIC) and 7 × 10⁻⁹ – 13 × 10⁻⁹ m s⁻¹ in nature (from Delaney et al., 1998). Note that as Martian values are less well known, the ranges of natural parameters defined here are for Earth only. The derivation and meaning of the Π-numbers are explained in Appendix A.

<table><thead><tr><th>Governing Parameter</th><th>Symbol</th><th>Model Min</th><th>Model Max</th><th>Nature Min</th><th>Nature Max</th><th>Unit</th><th>Dimension Function</th></tr></thead><tbody><tr><td>Volcano Height</td><td>H</td><td>0.01</td><td>0.05</td><td>1000</td><td>15000</td><td>m</td><td>[H] = L</td></tr><tr><td>Volcano Radius</td><td>R</td><td>0.017</td><td>0.087</td><td>2500</td><td>90000</td><td>m</td><td>[R] = L</td></tr><tr><td>Brittle Basement Thickness</td><td>B</td><td>0</td><td>0.15</td><td>0</td><td>40000</td><td>m</td><td>[B] = L</td></tr><tr><td>Brittle Material Density</td><td>ρ<sub>B</sub></td><td>1250</td><td>1300</td><td>1800</td><td>3200</td><td>kg.m<sup>-3</sup></td><td>[ρ<sub>B</sub>] = M.L<sup>-3</sup></td></tr><tr><td>Brittle Material Cohesion</td><td>C<sub>B</sub></td><td>12</td><td>120</td><td>1.00 × 10<sup>5</sup></td><td>1.00 × 10<sup>7</sup></td><td>kg.m<sup>-1</sup>.s<sup>-2</sup></td><td>[C<sub>B</sub>] = M.L<sup>-1</sup>.T<sup>-2</sup></td></tr><tr><td>Brittle Material Friction Coeff.</td><td>μ</td><td>0.75</td><td>0.95</td><td>0.1</td><td>1</td><td>–</td><td>–</td></tr><tr><td>Ductile Basement Thickness</td><td>D</td><td>0.005</td><td>0.8</td><td>20</td><td>270000</td><td>m</td><td>[D] = L</td></tr><tr><td>Ductile Basement Density</td><td>ρ<sub>D</sub></td><td>965</td><td>1000</td><td>1500</td><td>3500</td><td>kg.m<sup>-3</sup></td><td>[ρ<sub>D</sub>] = M.L<sup>-3</sup></td></tr><tr><td>Ductile Basement Viscosity</td><td>η</td><td>10000</td><td>40000</td><td>2 × 10<sup>13</sup></td><td>1 × 10<sup>20</sup></td><td>kg.m<sup>-1</sup>.s<sup>-1</sup></td><td>[η] = M.L<sup>-1</sup>.T<sup>-1</sup></td></tr><tr><td>Gravitational Acceleration</td><td>g</td><td>9.8</td><td>9.8</td><td>9.8</td><td>9.8</td><td>m.s<sup>-2</sup></td><td>[g] = L.T<sup>-2</sup></td></tr><tr><td>Timespan of Deformation</td><td>t</td><td>7200</td><td>18000</td><td>1 × 10<sup>12</sup></td><td>1 × 10<sup>14</sup></td><td>s</td><td>[t] = T</td></tr><tr><th>Pi-number</th><th>Symbol</th><th>Model Min</th><th>Model Max</th><th>Nature Min</th><th>Nature Max</th><th>Similarity</th><th>Comment</th></tr><tr><td>Volcano slope (R/H)</td><td>Π<sub>1</sub></td><td>1.73</td><td>1.74</td><td>1.25</td><td>12.70</td><td>Yes</td><td></td></tr><tr><td>Brittle basement thickness to Volcano Height (B/H)</td><td>Π<sub>2</sub></td><td>0.00</td><td>3.00</td><td>0.00</td><td>2.50</td><td>Yes</td><td></td></tr><tr><td>Ductile basement thickness to Volcano Height (D/H)</td><td>Π<sub>3</sub></td><td>0.10</td><td>9.00</td><td>0.01</td><td>33.75</td><td>Yes</td><td></td></tr><tr><td>Volcano–Basement density ratio (ρ<sub>B</sub>/ρ<sub>D</sub>)</td><td>Π<sub>4</sub></td><td>1.25</td><td>1.35</td><td>0.51</td><td>2.13</td><td>Yes</td><td>More commonly Pi#4 &lt; 1 in nature?</td></tr><tr><td>Friction coefficient (μ)</td><td>Π<sub>5</sub></td><td>0.75</td><td>0.95</td><td>0.1</td><td>1</td><td>Yes</td><td></td></tr><tr><td>Hubbert number (C<sub>B</sub>/ρ<sub>B</sub>.g.(H+B))</td><td>Π<sub>6</sub></td><td>0.004</td><td>0.98</td><td>0.0001</td><td>0.57</td><td>Yes</td><td></td></tr><tr><td>Ramberg number (ρ<sub>B</sub>.g.H.t/η)</td><td>Π<sub>7</sub></td><td>22.05</td><td>1834.56</td><td>0.18</td><td>5 × 10<sup>9</sup></td><td>Yes</td><td></td></tr><tr><td>Reynolds number (v.ρ<sub>B</sub>.D/η)</td><td>Π<sub>8</sub></td><td>7.0 × 10<sup>-12</sup></td><td>1.2 × 10<sup>-6</sup></td><td>3.8 × 10<sup>-26</sup></td><td>7.5 × 10<sup>-13</sup></td><td>No</td><td>Laminar regime in model and nature</td></tr></tbody></table>

For these sectioned models, *B* and *D* values are the average of ten thickness measurements distributed under the cone.

## 2.2. Digital Image Correlation

Overhead, time-lapse photographs of the deforming models were acquired at 1-minute intervals. We used a DSLR camera (18 megapixels, Canon EOS 550D) fitted with a 35.0 mm fixed lens and mounted ~1.3 m above the model basement. With a field of view of approximately 562 × 842 mm, the image resolution was approximately 3 pixels/mm.

To calculate horizontal displacements, velocities, and strains of the model surface, the time-lapse photographs were processed with the StrainMaster2D module of the commercially available digital image correlation (DIC) software DaVis by LaVision (Adam et al., 2005). This software uses a fast-Fourier transform approach to locate and track peaks in pixel amplitudes from one image to the next in the time-lapse series. DaVis also enables the use of an adaptive multi-pass correlation to improve the precision and spatial resolution of the displacement vector field. Here, we conducted a first iteration with 64 × 64-pixel windows and 75% spatial overlap, followed by two iterations with 32 × 32-pixel or 12 × 12-pixel windows and 50% overlap. To increase the signal-to-noise ratio, displacement vectors were stacked with respect to every fifth image of the sequence. Default cutoff values were used to filter and disable spurious vectors. Any subsequent gaps in the vector field were filled by bilinear interpolation. These steps improved the displacement vector precision to <0.1 pixels or <0.03 mm.

The DIC results are shown in map view in terms of cumulative horizontal displacement and strain (Holohan et al., 2013). The latter is shown as dilatation, i.e., surface area increase or decrease, which we interpret as overall horizontal extension (red) or contraction (blue), respectively, in the two-dimensional plane of view. For instantaneous cone emplacement, we show the horizontal velocity of the cone flank in terms of a representative velocity and a velocity time series. The representative velocity derives from the maximum horizontal displacement within a 40-minute window after emplacement of the cone. The velocity time-series comprise

values calculated from the maximum flank displacement in 10-minute windows after cone emplacement. For incremental cone emplacement, the velocity values derive from the maximum horizontal displacement within a 20-minute window after the start of each growth increment. The maximum displacement typically occurs near the base of the cone flank, and so it is from here that all velocity values are taken.

## 2.3. Scaling

We define eight dimensionless Π-numbers that describe both the system mechanics and the models' geometric, dynamic, and kinematic scaling to nature (Table 1). The dimensional analysis follows that of Merle and Borgia (1996) and De Matteo et al. (2022), albeit with modifications (see Appendix A).

Π₁, Π₂, and Π₃ represent the normalization of the cone radius, the brittle layer thickness, and the ductile layer thickness, respectively, to the cone height (a characteristic length). The numbers Π₄ and Π₅ represent two key parameters of the Mohr-Coulomb failure criterion: internal friction coefficient and cohesion, respectively. Π₆ represents the density contrast between the brittle and ductile layers, and hence the buoyancy forces. Π₇ accounts for the rate of loading from the growing cone to the rate of strain (i.e., flow) in the ductile layer (Ramberg number). Π₈ describes the flow regime (laminar or turbulent) within the ductile layer (Reynolds number). For the model to be physically similar to the natural system, each of these Π-numbers should have the same values in model and nature (Barenblatt, 2003), or at least be of the same order of magnitude (Merle, 2015), a condition that is generally satisfied (Table 1).

Given the dimensionless form of these parameters, the mechanical system explored here may be of relevance across a range of scales in nature (Merle, 2015). Its elements may thus have a variety of geological equivalents. The geometric scaling most closely approximates a stratovolcano overlying a weak or viscous (ductile) lithological lower unit (e.g., clay or salt) within a sedimentary basin (Borgia and van Wyk de Vries, 2003). As shown in Table 1, our models also have relevance for large shield volcanoes, however. On a basin scale, these may directly overly weak continental or ocean floor sediments (Poland et al., 2017). On a crustal scale,
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Structural outcomes and horizontal displacements in experiments of instantaneous cone growth. These representative results illustrate: (A) endmember spreading (model S-S63, *B*/*H* = 0.15, *D*/*H* = 0.1); (B) hybrid "spragging" (model S-S41, *B*/*H* = 0.1, *D*/*H* = 0.5); and (C) end-member sagging (model S-S55, *B*/*H* = 0.1, *D*/*H* = 0.7). The top row shows the original overhead photographs taken at the end of each experiment (c. 300 minutes after the start in each case). The bottom row shows fields of total horizontal displacement and related dilatational strain after 150, 70, and 20 minutes, respectively. The models shown here differ mainly in terms of their *D/H* ratios. *R/H* values are the same; *B/H* ratios are 0.10-0.15. (See Fig. 4A and Table S3 for exact parameters.) Dashed black circles in each image mark the initial basal diameter of the cone. Note that these are close-ups to see the structural detail; the diameter of the lateral model boundary is about 14–16 cm greater than the fields of view shown here.</figcaption></figure>

continental volcanoes may be linked to a mechanical system comprising a seismogenic (brittle) upper crust and an aseismic (ductile) lower crust (Watts and Burov, 2003). On a lithosphere scale, oceanic seamounts also overly a mechanical system comprising the seismogenic (brittle) upper part of an oceanic lithosphere, the aseismic (ductile) lower part of that lithosphere, and an underlying viscous asthenosphere (Watts and Zhong, 2000). The lower boundary of the model represents a transition to stronger and more rigid material, such as: crystalline continental basement or oceanic floor on a basin scale; the sub-continental mantle on a crustal scale; or the mesosphere on a lithosphere scale. For temporal scaling, a typical model duration of six hours corresponds to $3.4 \times 10^5$ to $3.4 \times 10^6$ years in nature (depending on the viscosity assumed for nature).

# 3. Results

## 3.1. Instantaneous cone growth

### 3.1.1. End-member and hybrid deformation styles

The model outcomes were classified into endmember and hybrid deformation styles based on the overall development of strain and displacement in each experiment. Fig. 3 shows three representative models with instantaneous cone growth; the evolution of each model is shown in the supplementary Figures S1–S3.

End-member spreading models (Fig. 3A) display the following structural features and displacement patterns: (1) subsidence of the central summit, with either minor horizontal contraction or, more typically, extension along polygonal fractures; (2) expansion

of the cone flanks, accommodated by normal and strike-slips faults in a typical pattern of leaf-like grabens and wedge-shaped horsts; and (3) contraction at or beyond the cone base. The basal contraction is commonly accommodated by a circumferential fold and thrust belt. Alternatively (Fig. 3A), it is accommodated by conjugate oblique or strike-slip faults. The horizontal displacements of the cone during spreading increase from the summit to the lower cone flanks and are consistently directed outward through time.

End-member sagging models (Fig. 3C) display: (1) Subsidence and contraction of the summit and cone flanks, accommodated by thrust faulting to produce flank terraces in a 'fish-scale pattern'; (2) a trough around the base of the cone; and (3) a peripheral ring of extension in the basement outboard of the cone accommodated by circumferential fissures, normal faults, and grabens. The horizontal displacements of the cone during sagging also increase from the summit to the lower cone flanks but are consistently directed inwards through time.

Hybrid deformation styles (Fig. 3B) span a range of outcomes between spreading and sagging endmembers. Hybrid 'spragging' cones display (to varying degrees): (1) contraction and vertical subsidence of the summit area of the cone, accommodated by terrace-defining thrust faults; (2) initial expansion and spreading of the cone flanks, accommodated by normal faults; (3) peripheral extension in the basement adjacent to the cone; (4) a radial expansion accommodated by radial fissures or grabens in the surrounding basement; and (5) late-stage sagging, as initial outward-directed displacement of the cone later becomes inward-directed (Figure S2).
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Model behaviors as plotted in dimensionless D/H vs. B/H space. (A) Data from instantaneously grown cones. Colored symbols represent data from this work; black filled symbols are data of Merle and Borgia (1996) (Experiments 1, 6, 7 & 8 in their Table 3; these have comparable initial parameters and boundary conditions to ours). Contours of average horizontal velocity on the cone flank were constructed through linear interpolation and by assuming that velocity values for D = 0 are zero (for data see Table S3). The green and yellow lines delimit field of spreading, sagging, and hybrid "spragging". Sagging occurs where $B/H > -0.18 * (D/H) + 0.3$, whereas spreading occurs where $B/H < -0.56 * (D/H) + 0.25$. (B) Results from incrementally grown cones. Here, the deformation behaviors during each increment of cone growth are plotted for differing growth paths. The B, D, and H values of the colored points are based on the (false) assumption of no deformation of the cone or basement and so are an approximation only. Black points with open symbols represent a test of the assumed vs. measured values; see the Results section. The velocity contours and field boundaries are for reference – they are the same as those in part (A).</figcaption></figure>

### 3.1.2. Control of geometric parameters

Fig. 4A displays the classification of instantaneously constructed models plotted in a phase diagram defined by the dimensionless ratios $\Pi_2 = D/H$ and $\Pi_3 = B/H$. Gravity-driven spreading results when both numbers are low, i.e., when the cone height is large with respect to the thicknesses of both the brittle layer and the ductile layer. Conversely, volcano sagging results when one or both numbers are high – that is, when the cone height is small with respect to either or both basement layer thicknesses. The domains of spreading or sagging are separated by a negatively sloped hybrid 'spragging' field. As discussed in Section 4, these new data are con-

sistent with, but also provide clearer context for, the past results of Merle and Borgia (1996) – see Fig. 4A.

#### 3.2. Incremental cone growth

##### 3.2.1. Evolution of deformation style during growth

As a volcano grows, and if the basement layer thicknesses do not increase, then both $B/H$ and $D/H$ should tend to zero with time. With a growing volcano, the system should thus move toward the origin of the $B/H-D/H$ parameter space, ideally along a linear path of constant $B/D$ ratio if deformation is homogeneous (Fig. 4B). As such the volcano deformation should pass consec-
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Deformation and displacement during incremental cone growth. (A) Initial endmember sagging; (B) subsequent hybrid "spragging" behavior; (C) late-stage endmember spreading. The representative model shown here (S-Sinc04) has an initial B/D ratio of ~0.3 (see Fig. 4B and Table S2). The top row shows the overhead photographs taken at the end of each growth increment. The bottom row shows fields of cumulative horizontal displacement and dilatational strain. Again, note that these are close-ups to see the structural detail; the diameter of the lateral model boundary is about 14–16 cm greater than the fields of view shown here.</figcaption></figure>

utively through sagging, 'spragging,' and spreading regimes. We tested this prediction for three incremental growth paths, each represented by a differing initial B/D ratio.

The results of one of these incremental growth models are shown in Fig. 5. This model has a mid-range growth path of B/D = 0.3 (Fig. 4B), and its structures and displacement patterns are like those observed on the other paths. Early cone growth results in a clear sagging style, with inward-directed horizontal displacements generating central contraction and peripheral extension (Fig. 5A). After the subsequent growth increment, however, a hybrid deformation style is seen. Horizontal displacement vectors are initially outward-directed (Fig. 5B), and subsequently rotate to become inward-directed over the 45-minute repose period. During later growth increments, spreading dominates. Outward-directed horizontal displacements dominate and are greatest on the cone itself, resulting in leaf-like radial valleys (grabens) and triangle-shaped upstanding wedges (horsts), with a fold and thrust belt in the surrounding basement (Fig. 5C).

Fig. 4B represents this transition during experimental cone growth on the *B/H* vs. *D/H* phase diagram. For purposes of plotting, we assume that values of *B* and *D* do not change substantially over time because they cannot be measured directly during the experiments. Under this assumption, the growth paths B/D = 0 and the B/D = 0.3 appear to transition to spreading slightly earlier than expected from the instantaneous experiment data. We assessed the error arising from this assumption by stopping three growth experiments of the same initial basement configuration (B/D = 0.5) after the third, sixth, or ninth growth increment, and then quickly cross-sectioning those models to measure *B*, *D*, and *H* at the end of those increments (Fig. 6). The measured B/D ratio for the cross-sectioned models is 0.47 on average and is thus

close to the assumed value of B/D = 0.50. The cross-sections show that both basement layers progressively thinned during growth, however, such that (as expected) the measured values of *B/H* and *D/H* lie closer to the origin than we had assumed (Fig. 4B). Moreover, this shift toward the origin increases slightly with cone growth. An origin-wards correction of similar magnitude applied to incremental-growth on paths B/D = 0 and B/D = 0.3 in Fig. 4B would make the transition to spreading there more compatible with the threshold values seen in the instantaneous-growth experiments in Fig. 4A.

### 3.3. Horizontal velocity of the cone flank

For instantaneously emplaced cones, horizontal velocity of the cone flank generally increases as both the *B/H* and *D/H* parameters decrease (Fig. 4A). For the limiting cases of *B/H* → ∞ (extremely thick brittle layer) and of *D/H* = 0 (no ductile layer), we expect no measurable velocity. *B/H* has a stronger effect than *D/H*, giving a closer vertical spacing of the velocity contours. The greatest horizontal velocities occur in the spreading field; the lowest occur in the sagging field.

Over time, velocities on the flanks of instantaneously emplaced cones decrease exponentially (Fig. 7A). This trend occurs for all deformation styles. Although we could not establish a relationship between *B/H* and the exponent, that exponent generally decreases with increased *D/H*. This finding means that, in general, horizontal velocities on a spreading cone (low *D/H*) remain higher for longer following emplacement compared with those on a sagging cone (high *D/H*).

For incrementally emplaced cones, on the other hand, the maximum horizontal flank velocity over the course of an experiment
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Digitized cross-sections through incrementally grown cones. Models were sectioned after (A) three, (B) six, and (C) nine growth increments. The initial B/D ratio (i.e., before emplacement of the first growth increment) was 0.5 in all cases. Deposits associated with each growth increment are numbered and visualized by alternating black and gray layers. Note how the basement is thinned below the deforming cone.</figcaption></figure>

generally increases with time (Fig. 7B). This behavior is consistent with the overall increase in velocity toward the origin of the B/H-D/H plot (Fig. 4B). For the growth paths B/D = 1.4 and B/D = 0.3, the velocity increase is near-linear. For growth path B/D = 0, the maximum horizontal velocity initially increases rapidly during the earliest growth increments before declining and then increasing again slowly during later growth increments. Velocities are generally lower along the growth paths of higher B/D ratios (Fig. 7B). This observation is consistent with the effect of the associated higher B/H ratios (Fig. 4B). Finally, horizontal velocity following a single growth increment declines exponentially with time (Fig. 7C) in a manner similar to behavior after the instantaneous emplacement of an entire cone.

# 4. Discussion

## 4.1. The role of basement rheology and geometry in volcano deformation

In common with many past studies, our main geometric and rheological assumption is that long-term, gravitationally driven deformation of a volcanic edifice is controlled by viscous flow within a single, laterally extensive basement layer. We have undertaken a systematic exploration of the geometric parameter space wherein the thickness (D) of that ductile-viscous basement layer and the thickness (B) of any overlying brittle-plastic layer are normalized to the volcano height (H). One new finding here concerns the geometric control on sagging and spreading (Fig. 4A). Previous postulations for a boundary between spreading and sagging domains included: (i) a positively sloped line at B/D ~ 2 (Merle and Borgia, 1996); (ii) a vertical line at D/H ~ 0.3 (van Wyk de Vries and

Matela, 1998); and (iii) a combination of positively-sloped and horizontal lines at B/D ~ 2 and B/D = 0.3 (Borgia et al., 2000). Our more extensive dataset shows that spreading and sagging domains are separated by a gradual, negatively-sloped zone of hybrid behavior in D/H-B/H space. Moreover, the B/D ratio is not as critical a control as previously thought, since we show that models with the same B/D ratio can either spread or sag.

The hybrid deformation style is also a new finding here. A distinctive characteristic of the hybrid deformation style is that the maximum cumulative horizontal displacement can occur in the basement rather than on the cone flank. This behavior is seen also in past numerical model simulations of volcano-basement interactions that were interpreted to result in 'basement extrusion' (van Wyk de Vries and Matela, 1998), and in a recent combined analog and numerical study of the geodetic behavior of Vesuvius (De Matteo et al., 2022). Here we show that the hybrid model behavior also features a transition in time from spreading-dominated to sagging-dominated behavior after cone emplacement (Figure S2). This temporal shift in cone deformation post-construction is compatible with faults on Vesuvius indicating that the volcano was spreading in the past (De Matteo et al., 2022) and with decadal cGPS data (especially horizontal components) indicating that Vesuvius is currently sagging (De Martino et al., 2021). Note that this spreading-to-sagging behavior after cone growth is the opposite of the sagging-to-spreading transition observed during cone growth (Fig. 6). Overall, 'spragging' behavior probably represents dynamic instability arising in the geometric transition zone between endmember behaviors.

A synthesis of literature data for volcanoes showing geological or geomorphological evidence of long-term, gravitationally driven deformation on Earth and Mars is shown in Fig. 8. In line with
<figure><img src="image_11.png" /><figcaption>(B) Incrementally-grown cones</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(C) Schematic evolution</figcaption></figure>

<figure><img src="image_13.png" /></figure>

**Fig. 7.** The evolution of flank velocity with time on model cones. (A) Horizontal velocity evolution of instantaneously grown cones. Velocity values were calculated from the maximum horizontal displacement observed during the previous 10 minutes (see Table S4). Shown are data from three representative experiments, each with a different initial D/H ratio. Flank velocity shows a power-law decay with time after cone emplacement, with an increasingly negative exponent as D/H ratio increases. (B) Horizontal velocity evolution of incrementally grown cones. The velocity values are calculated from the maximum horizontal displacement observed c. 20 minutes after the addition of each increment of cone growth (see Table S5). Shown are data from three experiments, each representative of a different initial B/D ratio (see Fig. 4B). Flank velocity generally increases with time as the cones grow, although this behavior depends too on the B/D ratio. (C) Schematic evolution of incremental horizontal flank velocity during incremental cone growth, illustrating exponential velocity decay within each growth increment and after cessation of growth.
<figure><img src="image_14.png" /><figcaption>**Fig. 8.** *D/H* vs. *B/H* values estimated from the literature for a range of natural volcanoes on Earth and Mars. Lower plot is a close-up of the lower left part of the upper plot (black box). Note that *B/H* ~0 for lizuna and Iwaki volcanoes, but these data are plotted slightly above the *D/H* axis for clarity regarding the error bars. See Table S6 for data sources, assumptions, and uncertainties.</figcaption></figure>

model results, volcanoes thought to be undergoing spreading plot in the lower left of *B/H*–*D/H* space, whereas those showing evidence of sagging (flexure) plot to the upper right. The natural basement parameters are subject to uncertainty; many are estimates based on previous authors’ analyses of the surrounding geological or geophysical data. The natural examples also span basin-and lithosphere-scales. Lithosphere-scale examples are not as well represented by the model scaling values. Nonetheless, the natural data are consistent conceptually regarding the interplay of *B/H* and *D/H* in determining whether a volcano spreads and/or sags. In particular, the occurrence of only sagging at the lithosphere scale can be explained by high *B/H* and *D/H* ratios from the large thicknesses of the lithosphere and asthenosphere relative to edifice height.

## 4.2. Untested parameters and model limitations

Other dimensionless parameters that may affect the model outcomes are cone slope ($\prod_1$), density ratio ($\prod_4$), and ductile layer viscosity ($\prod_5$). Our cone slope value of $R/H = 1.7$ is compatible with stratovolcano slopes but lies at the steeper end of the natural range ($R/H = 1.5–22.6$). In a test, a cone emplaced instantaneously with a gentler slope ($R/H = 2.75$, $B/H = 0.14$, and $D/H = 2.29$) gave a spragging outcome instead of a sagging outcome. This result agrees with a recent study (De Matteo et al., 2022) that used still gentler cone slopes ($R/H = 5.55$, $B/H = 0.64$,

and $D/H = 0.40$ - their model 8) and that also gave a spragging outcome instead of a sagging outcome. Moreover, past numerical models (van Wyk de Vries and Matela, 1998) that employed a gentler slope ($R/H = 3.0$, $B/H = 0.5$, and $D/H = 1.0$) resulted in spreading of the cone, rather than sagging (Fig. 8). We expect, therefore, that decreasing the cone slope should shift the transition zone upward and to the right within the $B/H – D/H$ space. A density ratio of $\prod_4 < 1$—that is, lower than in our experiments—is more probable in nature (Saballos et al., 2013; Zucca et al., 1982). As for decreased H, decreasing the cone density would decrease mass and could inhibit spreading. Decreased viscosity might promote sagging over spreading (van Wyk de Vries and Matela, 1998). Thus, decreased cone slope, cone density, and/or ductile layer viscosity could shift the transition zone position in $B/H–D/H$ space. Future work could test this proposition.

Gravitational instability in our models arises from a single widespread ductile or viscous layer of finite thickness within the basement. We thereby neglect, for simplicity, the superposition of deformation fields rooted at multiple ductile layers in the basement, although this is likely to be critical for certain volcanoes such as Hawaii and Olympus Mons (Byrne et al., 2013; McGovern and Solomon, 1993). Nonetheless, our results enable the effects of each layer to be assessed independently (Fig. 8). Alternatively, basement weakness may be envisaged as a frictional surface subject to high pore-fluid pressures (Dieterich, 1988; Ruch et al., 2010). Models with low $D/H$ values approximate the structural
outcomes of such a low friction control, though the dynamics may differ in detail (Rosenau et al., 2009). Low friction or ductile layers may not underlie the entire volcano in nature, and a volcanic edifice may be buttressed by pre-existing topographic features, cases which may yield more unidirectional or bidirectional flank motions (Cecchi et al., 2004; Merle and Borgia, 1996; Platz et al., 2011).

Volcano instability may also relate to factors other than those simulated (Poland et al., 2017). One is loading from dense intrusive complexes located within or beneath an edifice (Borgia et al., 1992; Delaney et al., 1990). This intrusion-related loading may interact with basement-induced deformation, potentially leading to positive feedback in the case of volcanic spreading (Le Corvec and Walter, 2009). Hydrothermal systems may create zones of weakened rock mass that facilitate flank creep, subsidence, and/or collapse (Cecchi et al., 2004). Finally, regional tectonic structures may guide magma intrusion and edifice construction, and they may influence structure and directionality of flank displacement (Grosse et al., 2020; Kervyn et al., 2014; Norini and Acocella, 2011). Simulation of these factors is beyond the scope of the present study.

### 4.3. Transition from sagging to spreading during cone growth

Our study demonstrates a sagging-to-spreading transition over time during incremental cone growth (Figs. 4B and 5). Such a transition was predicted by some studies with instantaneous cone emplacement (Borgia, 1994; van Wyk de Vries and Matela, 1998), but was not reported from previous modeling of cone growth (Cecchi et al., 2004; Delcamp et al., 2008; McGovern and Solomon, 1993). Here, we consider how flow within the ductile layer may explain the three-dimensional development of displacement and structures as a volcanic edifice grows incrementally (Fig. 9).

At the start of cone growth, a thick ductile basement layer relative to cone height enables a strong downward component of flow in response to the loading from the cone (van Wyk de Vries and Matela, 1998) (Fig. 9A). Subsidence with inward-directed horizontal surface motions of the cone leads to concentric and radial shortening of the cone flanks and to formation of flank terraces, reverse faults, and/or folds. Beyond a trough surrounding the cone, upward and outward flow in the ductile layer creates a peripheral bulge and a zone of extension with circumferential fissures and/or normal faults. As the cone grows, the downward component of flow becomes inhibited by the finite thickness of the ductile layer (van Wyk de Vries and Matela, 1998), and the flow is forced laterally (Fig. 9B). In the last increments of cone growth, flow below the cone is dominantly lateral (Fig. 9C). Subsidence of the cone is now accompanied by outward-directed horizontal surface motions, leading to concentric and radial lengthening of the cone flanks and to normal or oblique-normal faulting. The horizontal flow component decreases toward the edge of the cone, which leads to radial contraction and circumferential extension there, giving rise to a circumferential belt of reverse faults and/or folds that is intersected by near-radially orientated normal or strike-slip faults.

Evidence in nature for such a transition may include uplift, rotation, inversion, and cross-cutting of peripheral structures from the earlier sagging-dominant phases. At Concepción volcano (Fig. 1), field exposures reveal early normal faults interpreted as related to an early phase of radial peripheral extension induced by volcano sagging (Borgia and van Wyk de Vries, 2003). These faults have been progressively rotated and then uplifted within peripheral folds and monoclines, which were interpreted as having resulted from later peripheral radial contraction induced by volcano spreading. Our experiments indicate that this apparent transition from sagging to spreading occurred as *D/H* and *B/H* ratios decreased during Concepcion's long-term growth (Fig. 9).

For many volcanic seamounts such as Hawaii, the role of the basement in long-term, gravity-driven deformation is likely com-

plex. Two controlling layers of importance include: (1) the asthenosphere and the lower aseismic part of the lithosphere (Watts and Zhong, 2000); and (2) the weak oceanic sediments upon which the edifice is constructed (Dieterich, 1988). Our models suggest that these shallow and deep controls will be complementary in the early stages of volcano growth, as high *D/H* values promote sagging for both. In later growth stages, the shallow-level control from the oceanic sediment layer likely transitions into the spreading domain, as the associated *D/H* ratio decreases and thus gives rise to spreading-related rift zones. Indeed, seamounts of height greater than 2,000 m commonly exhibit stellate morphologies and rift zones (Vogt and Smoot, 1984; Mitchell, 2001).

Oceanic sediment thickness today is typically 10–1000 m (Straume et al., 2019). For *B/H* ~ 0, our experiments predict onset of spreading at *D/H* ~ 0.4. For *D* = 100–1000 m, a critical height above which spreading should occur is *H* = 250–2500 m. Although the upper end of this height range agrees with natural observations of height vs. rifting, our models predict spreading for some edifice heights lower than apparently observed. This discrepancy may have several explanations. Firstly, many non-stellate seamounts occur at mid-ocean ridges (Mitchell, 2001), where sediment thickness is typically only a few 10s of m. Here, the absolute thickness of ocean sediment may be too little to act as an effective décollement (McGovern et al., 2015). The absence of major volcanic rifts at the Galapagos Platform (Fig. 1) may relate to this factor (McGovern et al., 2015). Secondly, ocean sediment rheology may be non-Newtonian (e.g., Bingham), such that a certain shear stress (and thus volcano height) must be exceeded before flow can occur. Thirdly, seamounts of height <2,000 m may spread, but the structural expression of that deformation is not apparent until later in the volcano's growth (i.e., at heights >2,000 m) when the spreading-related velocities are sufficiently high. In our models, normal faults and grabens on the cone flank are usually not clearly visible until the latest growth increments (Fig. 5), at which point *D/H* ~ 0.1 (see cross-section test data in Fig. 4B). For *D* = 100–200 m in nature, this critical model height would correspond to a height of 1000–2000 m, i.e., much closer to the height values above which stellate forms are observed in nature. Overall, therefore, observations of seamount height, shape, and rift zone development are compatible with the general evolution of our models of cone growth influenced by basement-controlled deformation. Future work should further test this link and the influence of other factors as noted above.

### 4.4. Flank velocity during and after volcanic edifice construction

Experiments featuring instantaneous emplacement of the entire cone reveal mainly post-construction displacement and velocity (Fig. 7A). As in past works (Cecchi et al., 2004; Delcamp et al., 2008), our experiments indicate that basement-controlled volcano flank velocities undergo a long-term exponential decay after major eruptive phases or after extinction. Experiments with incremental emplacement of the cone show syn-construction displacement and velocity (Fig. 7B). These new experiments predict that long-term flank velocity generally increases in a linear to non-linear manner over multiple increments of cone growth. The horizontal velocity of flank motion during and after model cone construction thus shows distinct long-term evolutions, which are influenced by the basement configuration.

In general, the *B/H* ratio has the greatest impact on the flank's horizontal velocity. Instantaneously emplaced cones show consistently lower flank velocities at higher *B/H* ratios (Fig. 4). This behavior is replicated in the cone growth experiments, whereby flank velocities across all growth increments are generally lower for higher *B/D* growth paths — and thus higher *B/H* ratios (Fig. 7B). The role of *B/H* ratio is probably a function of three factors.
<figure><img src="image_15.png" /><figcaption>**Fig. 9.** Schematic illustration of the transition from volcano sagging to spreading during incremental cone growth. Shown are displacements, structures, and related kinematics for: (A) the early sagging phase; (B) hybrid 'spragging' phase; and (C) the late spreading phase.</figcaption></figure>
First, the brittle layer exhibits a pressure-dependent increase in its strength (Poppe et al., 2021), and so is stronger when thicker. Second, a relatively thick brittle layer has a greater effective flexural rigidity, which inhibits basement flexing and folding (Merle and Borgia, 1996). Third, the thicker brittle layer distributes the weight of the cone more widely across the ductile layer, reducing the overburden stress gradient and thus the strain rate.

The D/H ratio has a weaker and less consistent effect on flank horizontal velocity. On one hand, instantaneously emplaced cones show a lower horizontal velocity (Fig. 4) and a more rapid exponential decay of that velocity (Fig. 7A) at higher D/H ratios. This behavior probably reflects increased vertical displacement due to downward flow with greater ductile layer thickness (van Wyk de Vries and Matela, 1998). On the other hand, cone growth experiments with a B/D = 0 growth path are characterized by an initial velocity increase followed by a decrease, before a resumption of velocity increase across later growth increments (Fig. 7B). This trend could be a result of initial subsidence at high D/H ratio creating a buttressing effect against spreading that is only overcome in the later growth increments.

Our results provide context for the interpretation of medium- to long-term geodetic data at actively spreading or sagging volcanoes, although, given the long timescales of large volcano construction (10⁴–10⁶ years), it is unlikely that these model predictions can be tested entirely by geodetic means. Broadly, there are close similarities between the spatial increase of velocity downward along the model cone flank during spreading and geodetic observations of horizontal flank velocity distribution at Kilauea and Etna volcanoes (Poland et al., 2017). The decadal, near-linear spreading rate at Kilauea (Delaney et al., 1998) could be a tangent to our predicted curve of a linearly increasing spreading rate with time, but the natural time series are too short to be definitive. Factors other than those in our model, e.g., deformation of subvolcanic intrusions, may be more relevant on the decadal time scale. Exponential decay of spreading rate following growth increments (i.e., eruptive phases) is not evidenced in the Kilauea geodetic data (Delaney et al., 1998). In the limit, if growth increment volumes are infinitesimal with respect to the established cone volume and occur in rapid succession, then any exponential velocity decay will be suppressed and only a linear increase will be visible. Testing our model predictions probably requires other geophysical and geological approaches, such as field mapping, seismic imaging, deep drilling, and radiometric dating of volcanic and sedimentary growth sequences (Collier and Watts, 2001).

# 5. Summary and conclusions

1. Basement controls on volcano sagging and spreading are experimentally characterized in a dimensionless two-parameter space. The latter is defined from brittle and ductile basement layer thicknesses normalized to cone height. A newly constrained zone of hybrid deformation ('spragging') occurs between end-member sagging and spreading fields.

2. A transition from volcano sagging to volcano spreading is observed for incremental cone growth. As a cone grows, its height relative to the basement layer thicknesses increases, such that the system's geometry passes within the dimensionless parameter space from the sagging field through the transition zone and into the spreading field.

3. Distinct evolutions of the horizontal velocity of the cone flank with time are observed during and after cone construction. During cone construction, flank velocity decays exponentially after each growth increment but, overall, flank velocity increases linearly to non-linearly across multiple growth increments. After cone construction, flank velocities decay exponentially. The rates of flank

motion are controlled by the ratios of basement layer thicknesses to cone height.

4. A synthesis of literature data for volcanoes showing geological, geomorphological, and/or geophysical evidence of long-term, gravitationally driven deformation on Earth and Mars shows first order agreement with our experiment results. The observation of only sagging at lithosphere scale can be explained as a consequence of high B/H and D/H ratios; only at the crustal or basin scales are these ratios low enough to produce spreading. Furthermore, the experimentally predicted transition from sagging to spreading during cone growth is compatible with structural observations at Concepción stratovolcano and with a height-linked transition from circular to stellate shape (plus rift zones) at volcanic seamounts.

5. Our study therefore provides a refined conceptual basis for understanding long-term, gravity-driven deformation of volcanoes on a wide range of scales (basin, crust, lithosphere) on terrestrial planets. Effects of factors not accounted for here, such as major magmatic intrusions, regional tectonics, edifice buttressing, edifice elongation, etc., can be (re-)examined from this basis.

# CRediT authorship contribution statement

Eoghan P. Holohan: Conceptualization, Data curation, Methodology, Supervision, Writing – original draft. Sam Poppe: Conceptualization, Data curation, Investigation, Methodology, Writing – original draft. Audray Delcamp: Investigation, Methodology, Writing – review & editing. Paul K. Byrne: Writing – review & editing. Thomas R. Walter: Resources, Writing – review & editing. Benjamin van Wyk de Vries: Writing – review & editing. Matthieu Kervyn: Conceptualization, Resources, Supervision, Writing – review & editing.

# Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# Data availability

Data will be made available on request.

# Acknowledgements

This manuscript greatly benefited from the comments of four anonymous reviewers and editor C. Petrone. We acknowledge O. Galland and Beanhive for fruitful discussions, and M. Rosenau for access to the DaVis software. During this work EPH received financial support from University College Dublin, GFZ-Potsdam, the European Regional Development Fund project AGEO - Platform for Atlantic Geohazards (Grant code: EAPA_884/2018) and the Irish Centre for Research in Applied Geosciences (iCRAG - grant code: 13/RC/2092_P2). SP was supported by an Aspirant PhD grant from FWO-Flanders (2014-2019), and since 2021 by ULAM postdoctoral scholarship (PPN/ULM/2020/1/00067) from the Polish National Agency for Academic Exchange (NAWA) and project 2020/37/K/ST10/02447 from the Norwegian Financial Mechanism 2014-2021 (NCN Poland). AD acknowledges an FWO-Flanders Postdoc grant. BvWdF acknowledges funding from the UNESCO International Geosciences Programme (Grant no. 692: Geoheritage for Resilience). MK acknowledges funding from the Research Foundation - Flanders (FWO Credit 1505212N and FWO project G029820N).
## Appendix A. Dimensional analysis

The simplified system we consider is entirely mechanical (i.e., thermal effects are neglected), and is thus described by three fundamental units of length (meters), mass (kilograms), and time (seconds) (Barenblatt, 2003). The measured or output parameters of the system, the velocity (v) and strain (ε), are a function of 11 governing parameters: the radius (R) and height (H) of the cone; the thicknesses of the brittle (B) and ductile (D) layers; the friction coefficient (μ) and cohesion (C) of the brittle materials; the densities of the brittle (ρB) and ductile (ρD) materials; the dynamic viscosity (η) of the ductile layer; the gravitational acceleration (g); and the time span of deformation (t) (Fig. 2A). One governing parameter (μ) is dimensionless, whereas 10 are dimensional with dimension functions of the form:

$$[P] = L^x \cdot M^y \cdot T^z,$$

where [P] is the dimension function of a given parameter, P, and L, M and T are the fundamental units of length, mass, and time, respectively. The dimension functions describe the factor by which the numerical value of each parameter changes between different systems of units (e.g. SI vs. Imperial) (Table 1).

Three of these governing parameters have independent dimension functions. This means that the dimension function of each of the remaining governing parameters can be formed as the product of the independent parameters' dimension functions raised to powers (Barenblatt, 2003). For example, the dimension function of viscosity η is the product of the dimension functions of ρD, D, and t (Table 1) raised to powers as follows:

$$[\eta] = [\rho_D]^1 \cdot [D]^2 \cdot [t]^{-2} = M \cdot L^{-1} \cdot T^{-1}$$

Thus, the dimension function of viscosity can be said to be dependent on those of density, ductile layer thickness, and deformation timespan. In contrast, the dimension functions of either density, ductile layer thickness or deformation timespan cannot be written as the product of powers of the two other dimension functions, and so these dimension functions are independent of each other. Essentially, if we choose any three of the dimensional governing parameters such that one of the three fundamental units (length, mass, time) is present in one but not the others, then the dimension functions of this subset of the governing parameters will be independent of each other. Thus, the number of independent dimension functions will be equal to the number of fundamental units (Merle, 2015; Merle and Borgia, 1996).

The Buckingham Π theorem states that a model's physical similarity to nature is described by n-k dimensionless 'Π'-numbers, where n is the total number of governing parameters and k is the number of governing parameters with independent dimension functions (Barenblatt, 2003). For the system considered here, this is 11-3=8 Π-numbers. The Buckingham Π-theorem does not prescribe how the Π-numbers are defined, however. One must consider the physical meaning and importance of variables in the chosen Π-numbers.

For the geometric scaling, we normalize the cone radius, R, the brittle layer thickness, B, and the ductile layer thickness, D to the cone height, H (a characteristic length) to obtain the Π numbers:

$$\Pi_1 = \frac{R}{H}; \Pi_2 = \frac{B}{H}; \Pi_3 = \frac{D}{H}$$

For the dynamic scaling, we take the governing equation for brittle failure of the 90:10 sand/plaster mix (Poppe et al., 2021) and for rocks (Byerlee, 1978) to be the Mohr-Coulomb criterion:

$$\sigma_s = \mu \cdot \sigma_n + C,$$

where σs and σn are respectively the shear and normal stresses on the failure plane at failure, μ is the coefficient of internal friction, and C is the material cohesion. Normalizing each side of that equation by

$$\sigma_n \approx \rho \cdot g \cdot H$$

yields two dimensionless numbers that relate brittle material strength to the applied gravitational forces: the coefficient of internal friction and the Hubbert number (Merle, 2015):

$$\Pi_4 = \mu; \Pi_5 = \frac{C}{\rho \cdot g \cdot H}$$

Additionally, buoyancy forces arise from the density contrast between the overlying brittle domain and the underlying ductile domain, which can be described by the dimensionless number:

$$\Pi_6 = \frac{\rho_B}{\rho_D}$$

For kinematic scaling, a key consideration is the growth rate of the cone vs. the strain rate arising from flow in the ductile layer. The growth rate of the cone (dV/dt) is given by:

$$\frac{dV}{dt} = \frac{\pi ((R_2)^2 \cdot H_2 - (R_1)^2 \cdot H_1)}{3 \cdot dt}$$

where V is the cone volume, H1 and R1 are respectively the cone height and radius before a growth increment, with R2 and H2 being the height and radius after the growth increment. The strain rate in the ductile layer is given by the equation for Newtonian flow as:

$$\dot{\epsilon} = \frac{\sigma_d}{\mu}$$

where $\dot{\epsilon}$ is the strain rate, $\sigma_d$ is the differential stress (or shear stress) acting on the fluid, and μ is the dynamic viscosity of the fluid. With the assumption that $\sigma_d \sim \rho_B \cdot g \cdot H_2$, one can define a dimensionless number relating the strain rate, cone growth rate, and volume of the cone after the growth increment as:

$$\Pi_7 = \frac{\dot{\epsilon} \cdot V \cdot dt}{dV} = \frac{\rho_B \cdot g \cdot (H_2)^2 \cdot (R_2)^2 \cdot dt}{\eta \cdot ((R_2)^2 \cdot H_2 - (R_1)^2 \cdot H_1)}$$

For the entire growth history of the cone, $H_1 = 0$, $R_1 = 0$, and $dt = t$, such that $\Pi_7$ simplifies to a relation between stress, viscosity, and time known as the Ramberg number (Merle, 2015):

$$\Pi_7 = \frac{\rho_B \cdot g \cdot H \cdot t}{\eta}$$

Finally, the flow regime within the ductile layer is described by the Reynolds number:

$$\Pi_8 = \frac{\rho_D \cdot v \cdot D}{\eta} = \frac{\rho_D \cdot D^2}{\eta \cdot t}$$

The dimensional analysis above largely follows that of Merle and Borgia (1996) and De Matteo et al. (2022), but with a few adjustments. First, we instead define all geometric parameters in terms of a characteristic length, H, (Borgia et al., 2000), which largely controls the load of the cone on its substratum. In this way, the Π-numbers defined as B/D in those works is simply the ratio of our Π2 and Π3. As shown by our new data, B/D influences model velocity, rather than whether a cone spreads or sags. B/D also influences the spacing of fractures and folds on and around the cone (Merle and Borgia, 1996). Second, compared to Merle and Borgia (1996), we consider the coefficient of internal friction as a dimensionless number on its own, because this is a critical stress ratio and a parameter within the Mohr-Coulomb equation. Third,
we include an additional dimensionless number, our Π₅, to scale cohesion – the other key parameter of the Mohr-Coulomb failure criterion (Merle, 2015). In this way we avoid the slight redundancy of the Π₆ of Merle and Borgia (1996). Otherwise our Π-numbers are equivalent to those in the previous studies.

## Appendix B. Supplementary material

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2023.118012.

## References

Adam, J., Urai, J.L., Wieneke, B., Oncken, O., Pfeiffer, K., Kukowski, N., Lohrmann, J., Hoth, S., van der Zee, W., Schmatz, J., 2005. Shear localisation and strain distribution during tectonic faulting - new insights from granular-flow experiments and high-resolution optical image correlation techniques. J. Struct. Geol. 27, 283-301.

Barenblatt, G.I., 2003. Scaling. Cambridge University Press, Cambridge.

Borgia, A., 1994. Dynamic basis of volcanic spreading. J. Geophys. Res., Solid Earth 99, 17791-17804.

Borgia, A., Delaney, P.T., Denlinger, R.P., 2000. Spreading volcanoes. Annu. Rev. Earth Planet. Sci. 28, 539-570.

Borgia, A., Ferrari, L., Pasquare, G., 1992. Importance of gravitational spreading in the tectonic and volcanic evolution of Mount Etna. Nature 357, 231-235.

Borgia, A., van Wyk de Vries, B., 2003. The volcano-tectonic evolution of Concepcion, Nicaragua. Bull. Volcanol. 65, 248-266.

Brooks, B.A., Foster, J., Sandwell, D., Wolfe, C.J., Okubo, P., Poland, M., Myer, D., 2008. Magmatically triggered slow slip at Kilauea volcano, Hawaii. Science 321, 1177.

Byerlee, J., 1978. Friction of rocks. Pure Appl. Geophys. 116, 615-626.

Byrne, P.K., Holohan, E.P., Kervyn, M., van Wyk de Vries, B., Troll, V.R., Murray, J.B., 2013. A sagging-spreading continuum of large volcano structure. Geology 41, 339-342.

Cecchi, E., van Wyk de Vries, B., Lavest, J.M., 2004. Flank spreading and collapse of weak-cored volcanoes. Bull. Volcanol. 67, 72-91.

Chiocci, F.L., Coltelli, M., Bosman, A., Cavallaro, D., 2011. Continental margin large-scale instability controlling the flank sliding of Etna volcano. Earth Planet. Sci. Lett. 305, 57-64.

Collier, J.S., Watts, A.B., 2001. Lithospheric response to volcanic loading by the Canary Islands: constraints from seismic reflection data in their flexural moat. Geophys. J. Int. 147, 660-676.

De Martino, P., Dolce, M., Brandi, G., Scarpato, G., Tammaro, U., 2021. The ground deformation history of the Neapolitan volcanic area (Campi Flegrei Caldera, Somma-Vesuvius volcano, and Ischia Island) from 20 years of continuous GPS observations (2000-2019). Remote Sens. 13.

De Matteo, A., Massa, B., Castaldo, R., D'Auria, L., James, M.R., Lane, S.J., Pepe, S., Tizzani, P., 2022. An integrated modeling approach for analyzing the deformation style of active volcanoes: Somma-Vesuvius case study. J. Geophys. Res., Solid Earth 127.

Delaney, P.T., Denlinger, R.P., Lisowski, M., Miklius, A., Okubo, P.G., Okamura, A.T., Sako, M.K., 1998. Volcanic spreading at Kilauea, 1976-1996. J. Geophys. Res., Solid Earth 103, 18003-18023.

Delaney, P.T., Fiske, R.S., Miklius, A., Okamura, A.T., Sako, M.K., 1990. Deep magma body beneath the summit and rift zones of Kilauea volcano, Hawaii. Science 247, 1311-1316.

Delcamp, A., van Wyk de Vries, B., James, M.R., 2008. The influence of edifice slope and substrata on volcano spreading. J. Volcanol. Geotherm. Res. 177, 925-943.

Dieterich, J.H., 1988. Growth and persistence of Hawaiian volcanic rift zones. J. Geophys. Res., Solid Earth 93, 4258-4270.

Grosse, P., Poppe, S., Delcamp, A., van Wyk de Vries, B., Kervyn, M., 2020. Volcano growth versus deformation by strike-slip faults: morphometric characterization through analogue modelling. Tectonophysics 781.

Holohan, E.P., Walter, T.R., Schopfer, M.P.J., Walsh, J.J., van Wyk de Vries, B., Troll, V.R., 2013. Origins of oblique-slip faulting during caldera subsidence. J. Geophys. Res., Solid Earth 118, 1778-1794.

Kervyn, M., van Wyk de Vries, B., Walter, T.R., Njome, M.S., Suh, C.E., Ernst, G.G.J., 2014. Directional flank spreading at Mount Cameroon volcano: evidence from analogue modeling. J. Geophys. Res., Solid Earth 119, 7542-7563.

Le Corvec, N., Walter, T.R., 2009. Volcano spreading and fault interaction influenced by rift zone intrusions: insights from analogue experiments analyzed with digital image correlation technique. J. Volcanol. Geotherm. Res. 183, 170-182.

McGovern, P.J., 2007. Flexural stresses beneath Hawaii: implications for the October 15, 2006, earthquakes and magma ascent. Geophys. Res. Lett. 34.

McGovern, P.J., Grosfils, E.B., Galgana, G.A., Morgan, J.K., Rumpf, M.E., Smith, J.R., Zimbelman, J.R., 2015. Lithospheric flexure and volcano basal boundary conditions: keys to the structural evolution of large volcanic edifices on the terrestrial planets. Geol. Soc. Spec. Publ. 401, 219-237.

McGovern, P.J., Solomon, S.C., 1993. State of stress, faulting, and eruption characteristics of large volcanos on Mars. J. Geophys. Res., Planets 98, 23553-23579.

Merle, O., 2015. The scaling of experiments on volcanic systems. Front. Earth Sci. 3, 26.

Merle, O., Borgia, A., 1996. Scaled experiments of volcanic spreading. J. Geophys. Res., Solid Earth 101, 13805-13817.

Mitchell, N.C., 2001. Transition from circular to stellate forms of submarine volcanoes. J. Geophys. Res., Solid Earth 106, 1987-2003.

Moore, J., 1970. Relationship between subsidence and volcanic load, Hawaii. Bull. Volcanol. 34, 562-576.

Moore, J.G., Clague, D.A., Holcomb, R.T., Lipman, P.W., Normark, W.R., Torresan, M.E., 1989. Prodigious submarine landslides on the Hawaiian ridge. J. Geophys. Res., Solid Earth 94, 17465-17484.

Neri, M., Casu, F., Acocella, V., Solaro, G., Pepe, S., Berardino, P., Sansosti, E., Caltabiano, T., Lundgren, P., Lanari, R., 2009. Deformation and eruptions at Mt. Etna (Italy): a lesson from 15 years of observations. Geophys. Res. Lett. 36.

Norini, G., Acocella, V., 2011. Analogue modeling of flank instability at Mount Etna: understanding the driving factors. J. Geophys. Res., Solid Earth 116.

Platz, T., Munn, S., Walter, T.R., Procter, J.N., McGuire, P.C., Dumke, A., Neukum, G., 2011. Vertical and lateral collapse of Tharsis Tholus, Mars. Earth Planet. Sci. Lett. 305, 445-455.

Poland, M.P., Peltier, A., Bonforte, A., Puglisi, G., 2017. The spectrum of persistent volcanic flank instability: a review and proposed framework based on Kilauea, Piton de la Fournaise, and Etna. J. Volcanol. Geotherm. Res. 339, 63-80.

Poppe, S., Holohan, E.P., Rudolf, M., Rosenau, M., Galland, O., Delcamp, A., Kervyn, M., 2021. Mechanical properties of quartz sand and gypsum powder (plaster) mixtures: implications for laboratory model analogues for the Earth's upper crust. Tectonophysics 814.

Rosenau, M., Lohrmann, J., Oncken, O., 2009. Shocks in a box: an analogue model of subduction earthquake cycles with application to seismotectonic forearc evolution. J. Geophys. Res., Solid Earth 114.

Ruch, J., Acocella, V., Storti, F., Neri, M., Pepe, S., Solaro, G., Sansosti, E., 2010. Detachment depth revealed by rollover deformation: an integrated approach at Mount Etna. Geophys. Res. Lett. 37.

Saballos, J.A., Malservisi, R., Connor, C.B., La Femina, P.C., Wetmore, P., 2013. Gravity and geodesy of Concepcion volcano, Nicaragua. Geol. Soc. Am. Spec. Pap. 498, 77-88.

Straume, E.O., Gaina, C., Medvedev, S., Hochmuth, K., Gohl, K., Whittaker, J.M., Fat-tah, R.A., Doornenbal, J.C., Hopper, J.R., 2019. GlobSed: updated total sediment thickness in the World's Oceans. Geochem. Geophys. Geosys. 20, 1756-1772.

van Wyk de Vries, B., Matela, R., 1998. Styles of volcano-induced deformation: numerical models of substratum flexure, spreading and extrusion. J. Volcanol. Geotherm. Res. 81, 1-18.

Vogt, P.R., Smoot, N.C., 1984. The Geisha guyots - multibeam bathymetry and morphometric interpretation. J. Geophys. Res. 89, 1085-1107.

Walcott, R.I., 1970. Flexure of lithosphere at Hawaii. Tectonophysics 9, 435-446.

Watts, A.B., Burov, E.B., 2003. Lithospheric strength and its relationship to the elastic and seismogenic layer thickness. Earth Planet. Sci. Lett. 213, 113-131.

Watts, A.B., Peirce, C., Collier, J., Dalwood, R., Canales, J.P., Henstock, T.J., 1997. A seismic study of lithospheric flexure in the vicinity of Tenerife, Canary Islands. Earth Planet. Sci. Lett. 146, 431-447.

Watts, A.B., Zhong, S., 2000. Observations of flexure and the rheology of oceanic lithosphere. Geophys. J. Int. 142, 855-875.

Weijermars, R., 1986. Finite strain of laminar flows can be visualized in Sgm36-polymer. Naturwissenschaften 73, 33-34.

Zucca, J.J., Hill, D.P., Kovach, R.L., 1982. Crustal structure of Mauna Loa Volcano, Hawaii, from seismic refraction and gravity-data. Bull. Seismol. Soc. Am. 72, 1535-1550.