Earth and Planetary Science Letters

Many volcanoes exhibit evidence of long-term deformation driven by gravity. This deformation inﬂuences 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 (ﬂexure) 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 ﬁrst 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 ﬂank increases linearly or non-linearly during cone construction, but it decreases exponentially after construction ceases. Despite their geometric simpliﬁcation 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/).


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 defor-mation 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 4 -10 6 years), although in the short term they can manifest 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).
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 nearradially orientated faults with normal to oblique-normal or strikeslip 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, 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 basementcontrolled sagging to spreading can occur during volcano growth.

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 4 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 3 ) 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 3 at 45-minute intervals were made typically upon an initial seed cone of volume 3 cm 3 (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.

Table 1
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 −8 -1 × 10 −5 ms −1 in experiment (from DIC) and 7 × 10 −9 -13 × 10 −9 ms −1 in nature (from Delaney et al., 1998) For these sectioned models, B and D values are the average of ten thickness measurements distributed under the cone.

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 10minute 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.

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 andBorgia (1996) andDe Matteo et al. (2022), albeit with modifications (see Appendix A).
1 , 2 , and 3 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 4 and 5 represent two key parameters of the Mohr-Coulomb failure criterion: internal friction coefficient and cohesion, respectively.6 represents the density contrast between the brittle and ductile layers, and hence the buoyancy forces.7 accounts for the rate of loading from the growing cone to the rate of strain (i.e., flow) in the ductile layer (Ramberg number).8 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, 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.
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 × 10 5 to 3.4 × 10 6 years in nature (depending on the viscosity assumed for nature).

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 outwarddirected displacement of the cone later becomes inward-directed (Figure S2).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".

Control of geometric parameters
Fig. 4A displays the classification of instantaneously constructed models plotted in a phase diagram defined by the dimensionless ratios 2 = D/H and 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.

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-  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.
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.

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 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.
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.

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 transition in from spreading-dominated to sagging-dominated behavior after emplacement (Figure S2).This temporal shift in cone deformation post-construction is compatible with 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  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.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.

Untested parameters and model limitations
Other dimensionless parameters that may affect the model outcomes are cone slope ( 1 ), density ratio ( 4 ), and ductile layer viscosity ( 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 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.

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).
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.

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 longterm 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.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 4 -10 6 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).

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.

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.
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: 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: 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: 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 σ n ≈ ρ.g.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): 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: 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: where V is the cone volume, H 1 and R 1 are respectively the cone height and radius before a growth increment, with R 2 and H 2 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: where ε is the strain rate, σ d is the differential stress (or shear stress) acting on the fluid, and μ is the dynamic viscosity of the fluid.With the assumption that σ d ∼ ρ B .g.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: For the entire growth history of the cone, H 1 = 0, R 1 = 0, and dt = t, such that 7 simplifies to a relation between stress, viscosity, and time known as the Ramberg number (Merle, 2015): The dimensional analysis above largely follows that of Merle andBorgia (1996) andDe 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 5 , 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 6 of Merle and Borgia (1996).Otherwise our -numbers are equivalent to those in the previous studies.

Fig. 1 .
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).

Fig. 2 .
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).

Fig. 4 .
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 Table3; these have comparable initial parameters and boundary conditions to ours).
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).

Fig. 7 .
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 TableS4).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 TableS5).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.

Fig. 8 .
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 Iizuna 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.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 basinand 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.

Fig. 9 .
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.
. 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.