Lava Delta Formation: Mathematical Modeling and Laboratory Experiments

We analyze the dynamics of evolving lava‐fed deltas through the use of shallow‐layer mathematical models and analog laboratory experiments. Numerical and asymptotic solutions are calculated for the cases of planar and three‐dimensional flows fed by a point source upstream of the shoreline. We consider several modes of delta formation: a reduction in the driving buoyancy force; an enhanced viscosity of the submerged material; and the production of a granular subaqueous platform, over which a subaerial current can propagate. These modes of delta formation result in different behaviors. Under a steady supply of fluid upstream, the buoyancy‐driven case develops a solution with a steady subaerial delta and a subaqueous current which propagates at a constant speed, while the granular platform model extends the delta indefinitely. We determine a late‐time power‐law relation for the shoreline extent with time in this case. When the viscosity contrast is large, the model with an enhanced subaqueous viscosity is shown to mimic the initial dynamics of the granular platform model, but ultimately reaches a steady shoreline extent at sufficiently late times, as for the buoyancy‐driven model. The distinct behaviors of these models are further illustrated through laboratory experiments, utilizing the gelling reaction of sodium alginate solution in the presence of calcium ions as a novel analog for the abrupt rheological changes that occur when lava makes contact with water. These experiments provide quantitative verification of the buoyancy‐driven model in the absence of the reaction, and demonstrate the effects of a subaqueous platform qualitatively in its presence.


Introduction
When a volcanic lava flow enters an ocean or lake, it extends the pre-eruption shoreline (e.g., see Figure 1), forming a "lava delta" (Bosman et al., 2014;Jones & Nelson, 1970;Rodriguez-Gonzalez et al., 2022;Soule et al., 2021;Umino et al., 2006).A combination of cooling induced rheological changes and the reduction in gravitational driving forces control the morphology and evolution of the delta.Such events are common on volcanic islands on Earth and there is evidence that Olympus Mons, a large shield volcano located on Mars, was also formed as a volcanic island in an early ocean of liquid water (De Blasio, 2012;Hildenbrand et al., 2023).The relatively flat ground that a lava delta forms can be an attractive prospect for agricultural or infrastructural development.For example, the town of Tamaduste in El Hierro (Canary Islands) is built upon a lava delta formed during an historic eruption (Rodriguez-Gonzalez et al., 2022).However, such deltas pose several significant hazards during or shortly following an eruption.They are prone to collapse, which can in turn cause explosions due to the sudden quenching of molten rock and the production of steam and, in some cases, even tsunamis due to the displacement of large volumes of water by volcanic material (Poland & Orr, 2014).The identification in geological records of features associated with such water entry events has also been used to draw conclusions regarding the paleoenvironment, such as past sea levels and glacial extents (Jones & Nelson, 1970;Smellie et al., 2013).
When discussing the dynamics of lava delta formation, it is convenient to make distinctions between cases in which the lava flow becomes submerged as a coherent current and those in which the volcanic material is fragmented by hydrovolcanic processes (Moore et al., 1973;Soule et al., 2021).The former situation leads to a subaqueous gravity current after crossing the shoreline.In the latter, the fragmentation of the lava produces loose granular debris of variable grain size that fills in and builds up the seabed topography.The resulting granular platform then becomes over-ridden by a coherent subaerial lava flow "cap."The platform and cap are analogous to foreset and topset bed formations in river deltas and are sometimes referred to as such in the volcanological literature (e.g., Skilling, 2002;Smellie et al., 2013).In this paper we refer to the subaqueous foreset bed simply as a "platform," and the subaerial topset bed as a "cap," for conciseness and clarity.A common instability of the lava delta is thus associated with the presence of a heavy subaerial lava cap above a loosely packed, unconsolidated granular platform which can ultimately lead to its gravitationally driven collapse (Poland & Orr, 2014).The pair of images in Figure 1 show a delta built during a lava flow on Kilauea volcano (Hawai'i) in 2017, with plumes of steam erupting from where hot lava interacts with sea water, leading to an abrupt, steep edge to the subaerial lava flow.
In this paper we report mathematical models and an experimental study of the dynamics of lava delta formation in these different regimes.Our theoretical model is built upon a dynamical description that exploits the thinness of the flowing material relative to its streamwise extent, and balances the driving streamwise gravitational forces with flow resistance.Such shallow layer models are common in environmental applications of fluid dynamics; for example, see the review of Griffiths (2000) for a discussion of models of lava flows, while Schoof and Hewitt (2013) review shallow ice models of glacier flow.The modeling framework may also be extended to immiscible fluid flows within porous media (e.g., see Bear, 1972) and materials that exhibit a yield stress (e.g., see Balmforth et al., 2014).The commonality between these models and the theoretical developments presented below is that the flow is predominantly parallel to the underlying impermeable bed, so that fluid acceleration perpendicular to the underlying bed is negligible to leading order and there is a balance between the gradient of the pressure and the gravitational forces.Then, parallel with the bed there is a balance between streamwise pressure gradients and viscous processes.Our theoretical model, formulated in Section 2, follows conventional approaches for free-surface viscous gravity currents (Huppert, 1982).However, we modify the model to include either a change in buoyancy on submergence (Section 3), an order of magnitude increase in viscosity of the submerged fluid (Section 4.1), or the instantaneous creation of a granular platform at contact with water (Section 4.2).We integrate the theoretical models forwards in time, to reveal the development of the submerged currents or lava deltas in response to a sustained source upslope.We also explicitly construct steady states and establish selfsimilar regimes for the unsteady dynamics.Our experiments, described in Section 5, utilize an aqueous solution of sodium alginate, which undergoes a gelling reaction on contact with a bath of calcium chloride solution, as a novel analog for lava undergoing rapid rheological evolution.The experiments agree well with the theoretical model in the case of buoyancy driven inflation (i.e., in the absence of the gelling reaction) and illustrate the qualitative behavior of the platform-building model when the gelling reaction occurs.The majority of the paper considers two-dimensional currents, which are either planar or constrained in the cross-slope direction.We extend the model to explore the lateral spreading of an unconstrained, three-dimensional current fed by a point source of material upstream of the shoreline in Section 6, presenting numerical solutions, an analysis of the novel selfsimilar behavior that arises at late times, and complementary experiments.In Section 7 we conclude, arguing that the simple models developed in the paper provide some key insights into the dynamic evolution of lava-fed deltas, and provide some directions for further work.

Mathematical Formulation
We model the lava flow as a shallow, slow, gravity-driven current of viscous fluid on a slope of constant inclination, α.We employ Cartesian coordinates in which x measures the horizontal distance from the initial shoreline, z is the vertical distance above the water level, and y is measured across slope from some reference position (see Figure 2).The general form of the evolution equation for such a shallow-layer model, derived from conservation of mass and momentum and assuming a shallow slope angle, tan α ≪ 1, is given by (e.g., see Oron et al., 1997) where h is the thickness of the current, S is a flux factor that depends on the geometry and the material properties, e x is the unit vector in the x-direction, and ∇ H is the horizontal gradient operator ∇ H = (∂ x ,∂ y ).The ambient fluid (air or water) is treated as inviscid.The shallow slope approximation, tan α ≪ 1 corresponds to α ≪ 45°.Klein et al. (2023) carried out a broad study of the geomorphology of 47 volcanic islands and found ranges of 4.3°-13.1°and 2.6°-10.9°for the average onshore and offshore slopes, respectively.Thus this assumption is good for many volcanic settings, but is less appropriate for the steeper volcanic slopes.Even for 13.1°, we have tan α = 0.233 which represents a reasonably shallow slope.

Planar Geometry
For a planar current, fed from far upstream (x → ∞) by a line-source of constant flux, Q, per unit length in the ydirection, there is no y-dependence and the problem reduces to (2) For a Newtonian fluid of density ρ and viscosity μ we have (e.g., Huppert, 1982;Jeffreys, 1925), where g′ is the reduced gravity g ′ = g ρ ρ a )/ρ, with ρ the density of the viscous current and ρ a the density of the overlying ambient which changes from air to water when the current becomes submerged.The flux per unit width at uniform current thickness, Q = S tan α, resulting from Equations 2 and 3, corresponds to the Jeffreys equation (Jeffreys, 1925), Q = ρg′h 3 sin α/3μ, which gives the volume flux per unit width for a lava flow of width significantly greater than its depth (d ≫ h).These expressions are equivalent at leading order, since tan α ∼ sin α for shallow angles tan α ≪ 1.

Hele-Shaw Geometry
For viscous flow within a narrow channel (i.e., a Hele-Shaw cell) of width d, the flux factor becomes (e.g., Huppert & Woods, 1995).If the width d and typical flow depth remain much smaller than the characteristic length scale along the channel, one can generalize this flux factor so that S = S(h, d) (e.g., see Batchelor, 1967), with the factor reducing to the Hele-Shaw limit, Equation 4, when d ≪ h, and to the planar limit, Equation 3, if d ≫ h.This solution for a finite-width channel has been used in the context of lava flows which self-channelize through the formation of solidified levees (Lev & James, 2014;Sakimoto & Gregg, 2001;Tallarico & Dragoni, 1999).For the most part, however, we avoid this generalization and consider only the planar and Hele-Shaw limits.The planar regime (d ≫ h) is typical of many lava flows, even in the case that they self-channelize.For example, Chevrel et al. (2013) report the dimensions of an Icelandic lava flow, giving a width d ≈ 100 m, and a height h ≈ 5 m, in a channelized region of the flow.Similarly, Dietterich et al. (2021) report the dimensions of lava flows during the 2018 eruption of Kilauea, with a typical depth of 5 m, and widths ranging from 30 to 160 m.In particular, the narrow channel of the Hele-Shaw geometry is not of real relevance to volcanological lava flows, but is useful here to interpret convenient analog laboratory experiments.Although not considered in this paper, in the regime in which the viscous current is positively buoyant in the denser ambient, the finite-width and Hele-Shaw geometries are of relevance to laterally confined marine ice sheets, where lateral buttressing plays a significant role in the dynamics (Pegler et al., 2013).

Scaling and Non-Dimensionalization
The vertical length scale of the problem is set by the flow depth far upstream from the shoreline, H.This is related to the volume flux (per unit width), Q, and the inclination, α, via We may then define a horizontal length scale, L, and time scale, T , by After removing dimensions, the governing Equation 2 becomes ∂ ĥ where the hats identify dimensionless variables and Ŝ ĥ) ≡ S H ĥ) S(H) . (8) For the different cases specified above, and in the absence of rheological change, we have with n = 3 for a planar current and n = 1 in the Hele-Shaw cell.In Equation 9, the factor g ′ / g ′ ∞ incorporates any spatial variation in the reduced gravity from the value far upstream (identified by the ∞ subscript).Later we will also consider rheological change as the current enters the water, for which we will give the form of Ŝ explicitly.
We assume that the currents are initiated far upstream of the shoreline and initially propagate down the dry slope before first entering the water at x = 0.In terms of the dimensionless variables, we therefore impose ĥ → 1 as x → ∞.In the absence of any material changes as the current enters the water, the increased buoyancy of the viscous fluid in the denser ambient fluid (water compared to air) results in a reduction of the driving force and hence in the speed of the current.A natural consequence is that the current must inflate to conserve mass.This buoyancy-driven inflation corresponds to encoding the buoyancy change in the reduced gravity, g′, and is explored in Section 3.Alternatively, upon making contact with the water, the rapid cooling may prompt rheological changes and solidification to locally arrest flow, or, in combination with the generation of steam, may result in hydrovolcanic explosions that produce avalanching granular debris.In either case, a platform of almost stationary material can build up from the seabed to the water-level, on top of which the subaerial gravity current continues to propagate.We model both mechanisms of platform formation in Section 4.

Buoyancy-Driven Inflation
Since the density of air is much less than the density of the lava, we have g ′ ∞ = g and thus, in the regions where the current is submerged (i.e., where x ĥ > 0), the factor g ′ / g ′ ∞ is given by 1 R with R = ρ w /ρ and ρ w is the density of the denser ambient (i.e., water).We therefore set where Θ(z) is the Heaviside step function which takes the value 1 when z ≥ 0, and 0 otherwise.
The density of basaltic melt at atmospheric pressure is typically around 2.7 × 10 3 kg/m 3 (Lesher & Spera, 2015;Murase & McBirney, 1973;Stolper & Walker, 1980), which corresponds to a value of R around 0.4.However, molten lava is often vesicular, with gas bubbles resulting in a lower density.Close to the volcanic vent, the vesicularity can be greater than 70%, but this decreases rapidly with distance from the vent, with a vesicularity of 20%-40% reported at a distance of 10 km (Cashman et al., 1994;Dietterich et al., 2021).Thus R ≈ 0.5 for a vesicularity of 20%, but could even exceed 1 for the near-vent vesicularities exceeding 70%.The case in which a viscous current becomes positively buoyant in the ambient (R > 1) is relevant to the flow of marine ice sheets, and has been studied in some detail (Kowal et al., 2016;Pegler et al., 2013;Robison et al., 2010).In this paper we will assume the water entry is sufficiently far from the vent that the current is not buoyant in water and thus remains in contact with the bed throughout the motion.Thus we have a typical range 0.5 < R < 1.In the absence of a more constrained value, where a specific value is required, we either take R = 0.7 or the value appropriate to our experimental materials, R ≈ 0.83.
A numerical solution for the buoyancy-driven inflation of a planar viscous current (n = 3) with R = 0.7 is displayed in Figures 3a and 3b at various instances.The current propagates down the dry incline and, on reaching the shore line, inflates and continues its progress, approaching a steady state.The corresponding evolution of the front and shoreline positions is shown in Figure 3c (measuring t from the moment that the current reaches the waterline).In Figure 3d, a numerical solution is also shown for R = 0.7 in the Hele-Shaw geometry (n = 1) for comparison, showing an enhanced effect of the change in buoyancy, and a less curved nose to the current.Details of the numerical method used to integrate Equation 7 with Equation 10 are given in Supporting Information S1.

Steady-State Delta Lengths
In the steady state reached after buoyancy-driven inflation, the dimensionless thickness of the current tends to constant values, 1 and ĥ+ , upstream and downstream of the shoreline respectively.These are connected by conservation of mass, which demands that and dictates the dimensionless downstream height as a function of the buoyancy parameter, R. We then define a measure of the dimensionless length, L of the steady-state delta, relative to the shoreline of the current without any buoyancy inflation (which would be at a dimensionless length of unity), via This length turns out to correspond precisely to the position of the new waterline in view of the detailed depth profile of the steady solution (see Section 3.2).
In a planar (n = 3) or Hele-Shaw (n = 1) geometry, we thus have The delta length therefore becomes unbounded as R → 1, corresponding to the elimination of any density difference and the removal of the gravitational driving force (cf.Robison et al., 2010).

Steady-State Depth Profiles
The depth profile of the steady state follows from solving the constant flux condition, In each region, ĥ > x and ĥ < x, there is a uniform thickness solution, of dimensionless heights, 1 and ĥ+ > 1, as discussed above.However, ĥ must be continuous at the shoreline where ĥ = x, and Ŝ ĥ) is an increasing function of ĥ.A consequence of this is that any depth variation over the submerged current is ruled out (i.e., ĥ ≡ ĥ+ in x > ĥ).The steady state therefore thickens only on its approach to the shoreline, and maintains constant depth after submergence.Consequently, the shoreline corresponds precisely to x = ĥ+ , and the length of delta is L, as above.
The depth profile therefore follows from solving only ĥn (1 d ĥ in x < ĥ+ .The implicit solution is for n = 3, and x = ĥ + log ( ĥ 1 ĥ+ 1 ) when n = 1.These predictions for the steady state profiles are included in Figures 3b and 3d.

Highly Viscous Platform
If the lava does not fragment when placed in contact with water, but cools rapidly to significantly increase the viscosity, the local resistance to flow may create a high viscosity platform above which a lower viscosity, subaerial current propagates.To model this situation, we assume that contact with water instantaneously increases the viscosity of the fluid, leading to a two-layer structure to the current (see Figure 2b).The upper portion of the current, in z > 0, has the original viscosity; that below the water level (z < 0), has an enhanced viscosity given by μ/ɛ (ɛ < 1).This is similar to the buoyancy-driven case (Section 3), in that an increase to the viscous resistance of the submerged current decreases the mobility encoded through Ŝ ĥ) .However, the two-layer structure that arises in this scenario results in a different mathematical formulation.Following the usual analysis for slow viscous gravity currents, we may compute the velocity profile of the two-layer current in the planar problem in the region x > 0. Provided the current is not fully submerged (so h > x), the dimensional horizontal velocity is given by Journal of Geophysical Research: Earth Surface

10.1029/2023JF007505
If the current is fully submerged and h < x, on the other hand, the viscosity is enhanced throughout the depth, gravity is reduced by buoyancy, and Computing the flux and converting to dimensionless variables as before, we find the dimensionless evolution equation in Equation 7, but with a new flux function.Similarly, we may repeat the calculation with Hele-Shaw geometry.Both cases are captured in (1 R)ε ĥn when ĥ < x. (20) Here the Heaviside step function, Θ( x) , allows Ŝ to also capture the desired flux rule in the on-shore region ( x < 0) .We denote the shoreline position (where ĥ = x) by x f (t).The position at the front of the current, where ĥ = 0 and the current reaches the seabed, is denoted by xb (t).
Figure 4 shows a numerical solution for the planar (n = 3) viscous platform model with ɛ = 2 × 10 3 and R = 0.7, showing the initial development of a steep fronted viscous platform and overlying spreading viscous current, after reaching the shoreline.This regime of delta building is evidenced in the inset of Figure 4b and in Figure 4d by the close coincidence of the shoreline and front positions for a period of time after crossing the shoreline.As the weight of the delta builds, the enhanced viscosity of the platform eventually becomes insufficient to hold up the flow, and the platform slumps, becoming a subaqueous current that propagates with an enhanced thickness relative to the subaerial flow.In particular, as for the buoyancy-drive case, the late time behavior of this model is a steady-state shoreline position, with a subaerial traveling wave solution of constant velocity (see Figures 4b and  4c).This steady state is discussed in Section 4.1.1.
The thermal structure of lava flows can be complex.They typically have large Péclet numbers (measuring the significance of advective to conductive heat transport) (Griffiths, 2000;Griffiths & Fink, 1992a;Hyman et al., 2022), implying the development of thin surface boundary layers on cooling that can insulate molten lava in the interior.Nonetheless, flowing lava is often close to the solidus temperature, rendering the rheology sensitive to temperature changes (Griffiths, 2000;Hyman et al., 2022).Furthermore, cooling is known to be significantly enhanced for submarine compared to subaerial lavas.For example, Griffiths andFink (1992a, 1992b) evaluated the radiative and convective heat fluxes from lava flows in air and water, finding that the combined heat flux for submarine lava is approximately 20 times larger than for subaerial lava at the same surface temperature.Mitchell et al. (2008) suggest that cooling in shallow water could be enhanced even further by the occurrence of cracking of the rapidly cooling surface crust and the generation of steam, which does not occur at the higher pressures of the deep water flows considered by Griffiths andFink (1992a, 1992b).Under conditions of rapid cooling upon interaction with water, the glass transition is of particular relevance.The temperature dependence of lavas is then often parameterized according to the empirical Vogel-Fulcher-Tammann equation (Fulcher, 1925;Tammann & Hesse, 1926;Vogel, 1921): where A, B, and C, are fitting parameters that depend on composition (Giordano et al., 2008), and μ and T are measured in Pa s and K, respectively.For example, with parameters appropriate to basalt: A = 4.55, B = 5,000 K, and C = 600 K (e.g., see Giordano et al., 2008;Lev et al., 2012;Chevrel et al., 2018), and starting from a temperature of 1,000°C, we find a viscosity ratio of ɛ = 2 × 10 3 for a decrease of 180°C.Thus, in some situations the water-enhanced cooling of a lava flow may result in a bulk rheological change on relatively short timescales, which in the model above is treated as a significant, instantaneous viscosity increase at contact.In a more sophisticated model, a chilled crust could be incorporated into the model via a different functional form for the flux factor, S, and this functional dependence could evolve in response to an evolving temperature field.However, this modeling of the thermal evolution is beyond the scope of this paper and we instead focus on minimal models that exhibit delta building.

Steady States
As for buoyancy-driven inflation, the viscous platform model (with governing Equation 7, and flux factor Equation 20) admits a steady solution with a finite delta.Conservation of flux implies that, for the far-field subaqueous current height, ĥ+ , and delta length, L, we must have As discussed in Section 3.2, again the steady-state depth profile takes a form in which the submerged current is of uniform depth, ĥ = ĥ+ , and the current thickens upstream of the shoreline.An example of this predicted steady state solution in the submerged region of the current is included in Figure 4b.An immediate consequence is that the steady-state shoreline is again located at x = ĥ+ .Hence, for ɛ ≪ 1, the steady-state delta is very long (significantly larger than the corresponding delta built by buoyancy alone), and takes a significantly longer time to become established.In the limit ɛ → 0, the length of this steady-state delta diverges, corresponding to a situation in which the delta is built out indefinitely from the shoreline.As we will show in the following section, this behavior is also exhibited by a model in which a granular platform forms from volcaniclastic debris deposited by molten lava crossing the shoreline.

Granular Platform Model
When volcanic lava flows enter water, the generation of steam and rapid quenching can fragment the lava in a process known as hydrovolcanic, or hyaloclastic fragmentation (Smellie et al., 2013;Soule et al., 2021).We model this process as an instantaneous transformation of viscous molten lava into a granular medium, wherever the current comes into contact with the water.We further assume that the resulting granular material rapidly avalanches to rest down the submarine slope over a timescale that is significantly faster than that of the relatively slow viscous lava flow, producing a platform at a fixed angle of repose, θ r (see Figure 2c).
The angle of repose for such a granular platform could vary significantly.Froehlich (2011) measured the angle of repose for a range of coarse-grained rip-rap and reported values between 33°and 42°(excluding the results for smooth, round particles which do not well represent typical hydrovolcanic debris).It is known that the angle of repose can depend on the size of the grains (Elekes & Parteli, 2021), reduced gravity via the buoyancy effect of the water (Elekes & Parteli, 2021), and cohesion due to the presence of interstitial fluid (Samadani & Kudrolli, 2001).The effect of cohesion is likely to be small for the size of clasts forming the bulk of a typical lava delta bed (e.g., 5 50 cm reported by Soule et al. (2021)).In some cases, foreset bed slopes have been explicitly measured for lava deltas in the field.The studies of Moore et al. (1973), Tribble (1991), Soule et al. (2021), Di Traglia et al. (2018), and Pore ¸bski and Gradziński (1990), cover deltas on Hawai'i, Pico Island, Stromboli, and the South Shetland Islands.Their reported slopes mostly lie in the range 20°-40°.Thus a value of θ r taken in this range would be reasonable for a typical lava delta.
The packing fraction, ϕ, of the granular platform is also poorly constrained for lava delta beds.For hard, monodisperse spheres, random close packing is given by ϕ = 0.64.Loosely packed spheres will have a lower packing fraction than this, however hyaloclastic debris is polydisperse which will typically act to increase the packing fraction (e.g., Farr & Groot, 2009).The non-spherical nature of the grains will also affect the packing, and can generate both higher and lower packing fractions (Zou & Yu, 1996).Where we are modeling the platform as a granular medium, we therefore opt to use a packing fraction of 0.64, in the absence of a more accurate value.
Under the granular platform model, the evolution of the delta is then modeled as a viscous current propagating over a horizontal platform with a frontal boundary condition which couples the volume flux delivered to the advancing shoreline with the rate of advancement of the front.For the one-dimensional currents, if the front is at the position x = x f ( t ) , this boundary condition is given in dimensionless variables by In (d), the shoreline position, x f , is plotted logarithmically and over longer times; the red and blue dashed lines indicate the limiting behavior for early and late times, where is a constant that accounts for the slope of the underlying seabed, α, the angle of repose, θ r , and the packing fraction, ϕ, of the granular platform.This boundary condition amounts to the planar limit of the higherdimensional model of Appendix B, and corresponds to a volume constraint which dictates that the rate of increase of the volume of the platform is equal to the volume flux from the front of the current above (after correcting for porosity).Thus, if the underlying slope is significantly shallower than the angle of repose (tan α ≪ tan θ r ), then k reduces to the packing fraction, ϕ, to leading order.On the other hand, k is undefined for α ≥θ r , since it is impossible to build a static platform on a base which is inclined above the angle of repose.In the particular case when the platform has a steep slope (tan θ r ≫ tan α) and vanishing porosity, we have k = ϕ = 1.This case corresponds closely to a viscous platform model with ɛ ≪ 1, since in this regime the viscous platform is essentially static.Ultimately, however, for any ɛ > 0, the viscous platform eventually slumps under the weight of the material above it, and the models diverge.The equivalence of the models in an asymptotic sense, and the time scale on which they begin to diverge, are detailed in Appendix A.
The dimensionless evolution equation for the granular platform model becomes where the step function represents the transition from the sloping topography to the horizontal granular platform at x = 0. Supporting Information S1 again provides the details of the strategy we employ to numerically integrate Equation 25 with Equation 23.
Figure 5 shows a sample numerical solution of Equations 23 and 25 for the case of a planar current.Note that for this illustrative calculation, the angle of repose is taken to be much larger than the underlying slope (tan θ r ≫ tan α).Hence, the avalanched granular pile has a surface slope that becomes asymptotically large in the dimensionless scaling of the model, leading to a vertical line in the figure.The evolution of the front and shoreline positions appears in Figures 5c-5d.Notably, no steady state is reached and the shoreline grows indefinitely, with rate that gradually decreases due to the deepening of the platform.
Note that, because the current is generated by a constant flux, it propagates at constant speed over the dry incline before entering the water.In fact, the linear-in-time advance of the front continues for a short while beyond the moment that the current reaches the waterline.This continuation arises because x f ≪ 1, and so the front condition Equation 23 implies a negligible loss of material through the nose.The front is also relatively steep, with ∂ ĥ/∂ x ≫ 1, and so the change in slope is not felt to leading order.Consequently, the granular platform initially builds at the speed of the current on-shore, and x f ∼ t when t ≪ 1, as seen in Figure 5d.

Late-Time Behavior
At long times after the flow has entered the ocean and begun to form the lava delta, the motion enters a regime in which the front propagates as a power-law function of time.In this subsection we show how this phase of the motion is in self-similar form to leading order.
The volume flux condition at x = 0 requires ĥn+1 / x f ∼ 1, while the spreading at the front requires and ĥ = N where ζ = x/ x f .
On substituting this ansatz into the governing equation we find that where a prime denotes differentiation with respect to ζ.Thus, since n ≥ 1, the motion enters a quasi-static regime at long times, in which the time derivative is sub-dominant (i.e., the left-hand side of Equation 27 can be neglected).Integrating Equation 27 and applying the flux condition at ζ = 0 then yields In this quasi-static regime, the volume flux transported is constant along the current.Identifying the front as ζ = 1, and applying the spreading condition here, we find N = ̅̅ ̅ 2 √ and the solution and ĥ = ( This result for x f is compared with the planar (n = 3) numerical solution for k = 0.64 in Figure 5d and shown to capture the late time behavior accurately.

Experiments
Our experiments serve two purposes: first, we quantitatively validate the model for inflation due to a reduction in buoyancy.Second, we provide a qualitative experimental analog for the platform-building model.To achieve the latter we make use of an aqueous solution of sodium alginate and sugar, which reacts with the divalent cations in a calcium chloride solution to form a stiff hydrogel (Lee et al., 2017;MacKenzie et al., 2022).By introducing one solution into the other, we can therefore prompt an abrupt rheological change which mimics the transition that creates a platform in our model.

Experimental Setup
Figure 6 shows a sketch of the experiment: an alginate solution was fed into a 3 mm-wide rectangular tank with acrylic walls by a constant-flux pump.The channel was inclined at angle α to the horizontal, and the wall at the lower end retained a bath of either deionized (DI) water or a CaCl 2 solution.The former simulates buoyancydriven inflation because the current enters the bath and sinks to its base.In contrast, with calcium chloride solution, the current reacts on entering the bath to generate a gel that remained trapped at the shoreline by the side walls (we roughened the bottom and back wall of the channel with sandpaper in order to minimize any slip, although the front wall was left smooth for clarity of imaging).Freshly arriving unreacted alginate was then diverted over the gel to permit further contact with the bath and progressively build a shelf for the incoming current.Overflow from the end of the tank minimized changes to the water level in the bath as the alginate solution displaced the ambient fluid.
The experiments were imaged through the front wall of the tank using a camcorder positioned approximately 50 cm from the tank.The videos were recorded at a resolution of 1,920 × 1,080, with a field of view of approximately 45 cm × 25 cm, and a frame rate of 30 frames per second.The interface profile of the alginate current was extracted by thresholding according to the difference between the red and blue channels of the RGB image.The horizontal positions of the shoreline, x f (t), and the front of the current along the base of the tank, x b (t), were then determined and tracked as functions of time.

Materials
Per liter of the alginate and sugar solution, 8 g of medium viscosity alginic acid sodium salt (from MP Biomedical) was dissolved in 1 L of 1.20 g/ml sugar solution, producing an alginate and sugar solution of density 1.20 ± 0.01 g/ml.DI water was used in the preparation of the sugar solution to avoid any contamination by divalent cations in the alginate solution.In addition, a small amount (<0.1% wt/wt) of murexide indicator (ammonium purpurate, Sigma-Aldrich) was added to the solution, to dye it a reddish purple.This indicator undergoes a color change to a pale yellow hue in the presence of Ca 2+ ions, and has been used by Lee et al. (2017) to measure the degree of gelation of alginate.We therefore hoped to utilize this color change to indicate regions of reacted and unreacted alginate in our flow cell.However, we found that the rate of the color change reaction was too slow to track gelled regions in real-time during the experiment.Finally, a small quantity of commercially available glitter particles were added to the solution to act as tracer particles for visualizing the flow.The rheology of the alginate and sugar solution was tested in a Kinexus ultra + rotational rheometer (Malvern Instruments, Worcestershire, United Kingdom) using roughened, parallel plates of diameter 40 mm and gap 1 mm.The results are shown in Figure 7, and suggest that the rheology of the solution is weakly shear thinning and well fitted by the power-law model, τ = K γn , with n ≈ 0.882.Given the proximity of this shear index to the Newtonian value of n = 1, we choose not to include shear thinning in our theoretical model but note that this may be a source of some error in the comparisons with theory.Although not included in the theory of this paper, the weakly shear-thinning nature of the sodium alginate solution may strengthen its relevance as a lava analog, given that magmatic melts often exhibit weak shear thinning as a result of the crystal and bubble content of the suspension (e.g., Mader et al., 2013).For the unreacting experiments, the DI water had a density of 0.997 g/cm 3 .In the reacting experiments, the bath was produced by dissolving 20 g of calcium chloride dihydrate (Sigma-Aldrich) per liter of DI water, giving a density of 1.02 ± 0.01 g/cm 3 .Thus, R = 0.830 for the non-reacting experiments, whereas R = 0.850 for the reacting experiments.The solutions in both baths were dyed green with food coloring to provide additional contrast with the incoming alginate current.

Results
Figure 8 shows images from three experiments, indicating how the thresholding effectively identifies the outline of the alginate current (as well as the outline of an adjoining yellow measuring tape) and the qualitative differences between the unreacting and reacting cases.In the unreacting case (Figures 8a-8c) the current continues to propagate under the water, but thickens slightly after entering the water due to the buoyancy, as discussed in Section 3. In contrast, when the current enters the CaCl 2 bath (see Figures 8d-8f and 8g-8i), the alginate gels, extending the shoreline significantly through the progressive building of a shelf over which a subaerial current can propagate.The gel is just made visible by the subtle color change from reddish-pink to yellowish-orange associated with the murexide indicator.The shelf formed in this process is analogous to the viscous or granular beds formed in the theoretical model, in that it supports a roughly horizontal subaerial current or cap which provides the material required to extend the structure via volume flux over the evolving shoreline.Unlike in the theoretical models, however, where the granular or very viscous material fills the vertical column down to the sea floor (e.g., see Figure 5), the experimental shelf becomes locked in place by the sidewalls and creates an elevated platform for the overlying current.After an interval of "delta-building" in this manner, Figures 8d-8i illustrates two different types of failure of the shelf, as described later.The elevated shelf seen in these experiments is reminiscent of the experiments of Robison et al. (2010), Pegler et al. (2013), and Kowal et al. (2016), which model the dynamics of marine ice shelfs (including cases in which the shelf is laterally confined).In all of these cases an incident viscous gravity current forms a suspended shelf that extends over a bath of low viscosity fluid.In these previous experiments, however, the shelf was less dense than the ambient and is thus supported by buoyancy, whereas in the present work the shelf is denser than the ambient and is only elevated above the bath by the support of the side walls.These scales were then used to convert the experimental results into dimensionless units for comparison with the shallow-layer theory.The comparison in Figure 9 indicates that the theory is successful in reproducing the experiment.Some discrepancy can be observed, which could arise for a number of reasons.First, the channel in the experiment was not particularly narrow in comparison to the flow depth, suggesting that the Hele-Shaw flux limit taken in the shallow-layer theory was not accurate.Nevertheless, when the flux law was replaced by one applying for finite width (S = S(h, d), as mentioned in Section 2.1.2),and the solutions recomputed, the theoretical time series did not change significantly.Second, a Newtonian flux rule was used, even though the alginate-sugar solution is slightly shear thinning.Finally, some discrepancy arises when the current initially submerges ( t = 0) .The dynamics of this initial contact with the water could be being complicated by surface tension, non-shallow effects, and the motion of the ambient fluid, which are all ignored in the model.Given these uncertainties, the agreement in Figure 9 is perhaps surprisingly good.
Figure 10 shows the corresponding time series for the reacting experiments in Figures 6c-6i (panels b, c), together with two other tests with similar volume fluxes (panels a, d).Although the experiments were intended to be similar, the detailed dynamics turned out differently each time, due to the failures of the shelf occurring at different moments.In all four cases, the gelled shelf created by the reaction supports the current off the base of the tank, bringing the bottom of the current to rest while the shoreline continues to move outward.The shoreline evolution that results is qualitatively similar to that captured by our platform models.To emphasize this point, in Figure 10b, we have included a solution of the granular platform model using the Hele-Shaw flux rule and a fitted value of k = 0.1.The shoreline position of this solution matches that observed in the experiment until the gel shelf fails close to the original shoreline.We note that, although the low packing fraction implied by k = 0.1 is not appropriate to lava delta platforms, it does represent a reasonable value for these reacting experiments.This is because, in contrast to a granular platform, the gel shelf occupies a small volume below the water line which significantly reduces the effective packing fraction of the region below the cap.Alternatively, one can see this by considering the gel shelf to be a platform of unit packing fraction, ϕ = 1, but with an angle of repose close to 180°.In this case, the denominator of the geometric factor in the definition of k (24) becomes large.
Figure 10 also highlights how failures of the gel shelf invariably interrupt the delta-building phase before the shelf becomes very long.Moreover, such failures take one of two possible forms.In the first of these, a break occurs in the shelf close to its base at the bottom of the tank, due to the growing pressure of the fluid above (such events are marked by the letter A in the figure).After this failure, the incoming fluid is largely driven through the gap in the shelf, feeding a curved submarine current (see Figures 8e and 8f) and slowing the advance of the shoreline.We attribute the curving of the interface of the submarine current to the additional resistance provided by a thin layer of gel created as fresh alginate solution comes into contact with the CaCl 2 bath.
The second failure mode of the shelf, indicated by points marked B in Figure 10, occurs at the tip of the shelf.In this case, failure arises when new material arriving at the tip abruptly sinks into the bath rather than building the shelf out further.The shoreline then temporarily retracts and, in the case of panel (c), creates a sudden jump in xb if the sinking shelf material reaches the base of the tank in advance of the propagating submarine current.What prompts the sinking of fluid at the tip of the shelf, and this second mode of failure, is not clear.A slight widening of the channel or spatial variations in the strength of the formed gel could be responsible.Indeed the locations and times at which failure occurred varied significantly between the four examples.Typically, after this point, material sinking from the shelf merges with the submarine current and for the remainder of the experiment the fluid propagates down the slope as a single coherent current.
Although both modes of failure terminate the building of a floating shelf, the shoreline continues to advance at later times due to the thickening of the curved submarine current.This advance is slower, and more comparable to that driven only by buoyancy differences (see Figure 9).For the latter, we anticipate that a steady-state shoreline position is eventually reached (as predicted theoretically).The reacting flow may also converge to such a steady state after failure of the gel shelf, if the layer of gel at the interface arrests motion there and isolates the fresh alginate solution flowing in the current underneath.

Three-Dimensional, Unconstrained Currents
We now consider unconstrained flows fed by a localized point source on the slope.The current thus develops a non-trivial three-dimensional shape before and after intersecting the water.

Theoretical Model
A granular platform model for the three-dimensional case can be derived in the same manner as Section 4. The dimensionless governing equation is ∂ ĥ where the flux factor Ŝ = ĥ3 , the Heaviside step-function Θ( x) represents the transition from the sloping topography to the horizontal granular platform at x = 0, and Q is a flux source term, taken to be of the form representing a source with a Gaussian profile that delivers a unit volume flux concentrated over a region of radius O(δ) centered at the point x 0 (thus reducing to a point source of unit flux in the limit δ → 0).The shoreline is now described by an evolving curve on the x ŷ plane, corresponding to a horizontal contour of the free surface of the current (see Figure 2d).The volume flux through this contour must be normal to the shoreline, and therefore builds up an underlying granular platform in that direction.The progression of the shoreline is dictated by the solid fraction, ϕ, the angle of repose, θ r , the distance from the original shoreline, x, and the normal direction to the evolving shoreline, n = (cos β, sin β).The derivation of the boundary condition, based on a Coulomb model for the granular platform, is given in Appendix B. For a given point on the shoreline xs , ŷs ) , the front condition is where k is given by Equation 24, as for the planar case.The first factor on the left hand side of Equation 33 arises due to the dependence of the granular platform's geometry on the normal direction of the shoreline, with the implication that the rate of shoreline advance is slightly enhanced in the cross-slope direction, compared to the down-slope direction.In the regime in which the angle of repose is significantly larger than the slope, we have tan α ≪ tan θ r and the front condition reduces to which is essentially equivalent to the planar model, but with the shoreline advancing in a direction normal to the boundary, rather than simply in the x direction.
In principle, and in conjunction with the evolution Equation 31, the free boundary condition Equation 33 could be integrated into a numerical simulation of an unconstrained three-dimensional current.However, this implementation is practically challenging, requiring an adaptive domain governed by Equation 33, and beyond the scope of this paper.Instead, for numerical simulation we resort to the three-dimensional generalization of the viscous platform model, which is more straightforward to implement and becomes equivalent to a granular platform model in certain limits (namely tan α cot θ r → 0, ϕ → 1 and ɛ → 0, see Appendix A).In particular, when ɛ is small but finite, we have shown that the two models agree at early times but eventually diverge due to the slumping of the viscous platform.In Section 6.2 we derive a late time similarity solution for the granular platform model, and thus provide a description of the granular platform model at times beyond the validity of the viscous platform approximation below.
In dimensionless variables, the three-dimensional viscous platform model is where the viscosity ratio is again set by ɛ, and the flux function is Journal of Geophysical Research: Earth Surface In our simulations, the source term is given by Equation 32 with the parameters δ = 0.25 and x 0 = ( 4, 0).We initialize the computation with a pre-wetted film of thickness ĥ = h ϵ = 0.01.While this film is not especially thin, the flux within it is order 10 6 , much smaller than the imposed source flux and implying only a small contribution to the dynamics.For boundary conditions, we take ŷ = 0 to be a line of symmetry (allowing us to compute a solution for only half of the domain, ŷ > 0), set ĥ = h ϵ at the up-and cross-slope boundaries, and a free flux outflow condition, ∂ x ĥ = 0, at the down-slope boundary.All these lateral boundaries are chosen to lie at sufficient distances from the source that they play no role in the dynamics.The evolution equation is solved using the finite element method as implemented in FEniCS (Alnaes et al., 2015) and additional details of the computational method are provided in Supporting Information S1.
Two sample solutions are shown in Figure 11 for buoyancy-driven inflation (ɛ = 1; panels (a, b)), and a very viscous platform with ɛ = 10 3 (panels (c, d)).In both cases, we set R = 0.833, guided by experimental conditions.These two solutions indicate how viscous-platform building again produces a significantly larger delta than purely buoyancy-driven inflation.The viscous-platform delta also extends prominently across the slope because the material passing through the shoreline must displace the entire vertical water column; the deeper the water, the slower the advance.Vertical cross sections through the two solutions along the center line, ŷ = 0, are shown in panel (e), demonstrating qualitative similarity to the structure of the planar currents, as seen in Figures 3 and 4.

Late Time Similarity Solution for a Three-Dimensional Granular Platform
We can gain analytical insight into the shape and evolution of the three-dimensional delta via the construction of a similarity solution at late times ( t ≫ 1) .The dimensionless model equation in x > 0 is given by ∂ ĥ subject to a unit volume flux delivered to the delta from upstream of the original shoreline, x = 0.This demands that ĥ3 ∂ ĥ where y f ( t ) is the lateral extent of the delta at x = 0 and the delta is assumed symmetric about ŷ = 0.The advancing front of the delta is located at x = x s ( t) = xs ( t), ŷs ( t)) , and the conditions at the front are given by ĥ xs , ŷs where n = (n x , n y ) is a vector normal to the boundary and cos β = n x /‖n‖.
First we establish the gearing between the spatial coordinates ( x, ŷ) and time, t, in this regime.If the leading order dynamics are quasi-static, such that |∂ ĥ/∂ t|≪|∇ H ĥ3 ∇ H ĥ) ⃒ ⃒ , then we require that x ∼ ŷ.(The validity of this quasistatic regime is verified below.)The inflow boundary condition at x = 0, Equation 38, requires that ĥ4 ŷ/ x ∼ 1.
Finally, the frontal condition, Equation 39b, requires that k x2 / t ∼ ĥ4 / x.Thus, combining these scalings we deduce that ĥ ∼ 1, while x ∼ ŷ ∼ ( t/k) 1/ 3 .We seek a similarity solution of the form Y, and ĥ = H(X,Y), (40) and the width of the delta at x = 0 is given by y f = y 0 ( t/ k) 1/ 3 (thus Y = 1 at X = 0 in similarity variables).The constant, y 0 is at this stage undetermined, and is fixed by enforcing the boundary conditions.Immediately we can deduce the motion is quasi-static to leading order when tk 2 ≫ 1; in other words, this regime emerges at late times.
It is convenient to write the location of the advancing boundary of the delta in the self-similar regime as X = X s (Y) (which is possible since X ′ s ≤ 0); in this case, a normal vector is n = 1, X ′ s ) .Furthermore, we write U = H 4 /4 and note that the quasi-static problem is linear in U.The self-similar problem is therefore Laplace's equation The inflow arises from the viscously controlled flow down the plane, and at late times the width of this sustained current becomes steady, which means that it is O t 1/ 3 ) in terms of the self-similar variables for the delta Equation 40.Therefore at late times the inflow looks like a point source of unit flux at the origin, which we represent as a Gaussian profile of standard deviation Y i ≪ 1, so that The boundary conditions at the edge of the advancing delta in the self-similar regime are given by where r = tan α cot θ r .Finally, symmetry of the delta requires The numerical method used to obtain the solution to this similarity problem is given in Supporting Information S1.We use this method to calculate solutions for r = 0, 0.1 & 0.5.The case r = 0 corresponds to the simplified limit in which the angle of repose is vertical (or equivalently the basal slope is very shallow, α ≈ 0°).For an angle of repose of θ r = 40°, r = 0.1 corresponds to α ≈ 5°and r = 0.5 corresponds to α ≈ 23°(both to the nearest degree).At a shallower angle of repose, θ r = 20°, r = 0.1 corresponds to α ≈ 2°and r = 0.5 corresponds to α ≈ 10°.We find that y 0 = 1.746, 1.760, 1.856 (for r = 0, 0.1, 0.5, respectively) and we plot the shapes of the advancing boundaries in Figure 12a.Notably, aside from scaling via k, the similarity solutions do not differ significantly from the simplified setting in which the angle of repose is vertical and r = 0.The small increase in y 0 with r corresponds to the delta extending slightly further across-slope, due to the directional dependence of the granular platform geometry.
To compare the similarity solution with the full numerical solution for the viscous platform model, we first set r = 0 and ϕ = 1.Secondly, since the similarity solution is unchanged under a translation of the time variable, it requires some care to select the similarity solution time, ts , for a given time of the numerical solution, t.With this in mind, rather than matching the times directly, we choose instead to select the similarity time such that the volume in the delta is equal to that found in the numerical solution.Specifically, for the similarity solution, the total volume in the delta at ts is given by where the first term is the leading order contribution from the subaqueous granular platform, and the second term is the volume in the subaerial capping flow.The constant V a is determined by integration of the current height, H = (4U) 1/4 , over the subaerial domain, For ϕ = 1 and r = 0 this evaluates to V a = 1.86.For the numerical solution (and, for an assumed constant source flux, in the field) the total volume in the delta is given approximately by V = t t0 , where t0 is the time at which the flow first crosses the shoreline.Thus, to compare solutions, we determine the similarity time ts by solving ts + 1.86 t2/ 3 s = t t0 .Following this convention we compare the similarity shape of the delta with the full numerical solutions in Figures 12b-12d.For this simulation, the current was found to cross the shoreline at t0 = 8 and the values of ts were calculated correct to 2 significant figures.The similarity shape captures the full solution well until later times (panel (d)), when the slumping of the viscous platform results in divergence between the viscous and granular platform models, as discussed for the planar case in Appendix A.
It is potentially challenging to compare the geometry shown in Figure 12 to the deltas observed in volcanological settings, because the latter are subject to various additional factors including erosion by wave action and collapse of the unstable lava bench (Mattox & Mangan, 1997;Moore et al., 1973;Rodriguez-Gonzalez et al., 2022).However, at the very least, the aspect ratio (of lateral to downslope extent) predicted from the similarity solution is consistent with field observations which note that lava deltas can extend significant distances along the shoreline, and become wider than the typical width of the feeding onshore flow (Moore, 1975;Moore et al., 1973;Umino et al., 2006).For example, Moore et al. (1973) reported a lava delta formed on Hawai'i in 1971 with a lateral dimension of 1.5 km and a seaward extent of 450 m, which closely matches the aspect ratio of the similarity solution in Figure 12 (and given by 1/(2y 0 ) = 0.286).Similarly, the shoreline geometries reported by Umino  (2006) during the formation of a Hawaiian lava delta in 1990 are reasonably consistent with the prediction of the similarity solution.In addition to this predicted geometry, the similarity solution also provides the rate of advance of the shoreline.

Unconstrained Experiments
Our experiments for unconstrained currents were conducted with a tilted tank as sketched in Figure 6.The tank, however, was much wider than that used in Section 5, with a width of approximately 20 cm.The pipe from the pump was also positioned a small distance above the slope, so that the fluid was delivered vertically to approximate a localized source of flux.The camcorder was used to record a top-down view; a second camera took some additional images from the side.
Figure 13 shows snapshots of an unreacting experiment (panels (a)-(c)) and a reacting experiment (panels (d)-(i)).For the unreacting experiment, the current quickly submerges and propagates along the bottom of the tank, whilst the shoreline converges to steady state.By contrast, for the reacting experiment, the shoreline continues to extend with the front of the current, proceeding significantly beyond the final position reached in the unreacting case.The reacting current also loses symmetry, as seen in panels (f)-(i).In the absence of the constraining side walls, the "delta" built in this case is not due to the formation of an underlying elevated shelf, but rather by a gel sack that dams a pool of unreacted alginate.As for the shelf, however, delta-building again terminated when the sack failed, either by over-topping or by breaking near the base.In the example shown in Figure 13, this failure resulted in a submarine current which propagated continuously despite the ongoing formation of gelled material.The growth and release of the pool of unreacted alginate can be seen in the side views of the experiment shown in Figures 13j-13l.
These unconstrained experiments, therefore, provide the same analogy to the distinct modes of delta formation explored in the theoretical models and constrained experiments.Namely, the delta formed by buoyancy alone quickly attains a finite steady shoreline position, and the experiments featuring the gelling reaction exhibit a much more significant extension to the shoreline.We note that these experiments provide only a qualitative analogy for the buoyancy-driven and platform-building modes of delta formation, since the physical process involved in forming a gel sack does not closely represent delta formation on the scale of a volcanic lava flow.Nonetheless, the reacting experiments demonstrate interesting morphologies that mimic some smaller scale features of lava flows.The stalling of the flow front by the thin gel sack, followed by growth of the pooling material, and eventual breakout by rupture, is reminiscent of the pahoehoe style emplacement of lava via the inflation and rupture of individual flow toes (Hon et al., 1994;Rowland & Walker, 1990).Similarly, subaqueous lava often forms pillowy structures in which molten lava is insulated from the surrounding water by a solidified crust (Jones, 1968;Moore & Lockwood, 1978;Tribble, 1991).Growth of these pillows is reported to result from the formation of cracks which spread, exposing molten lava to the water and forming new crust on the exposed surface.

Conclusions
Motivated by the geophysical phenomenon of lava delta formation, we have explored theoretically the timedependent dynamics of a viscous gravity current entering an ambient water bath.First, we considered the case in which no rheological change occurs, but the current slows and inflates in response to a weakening of the driving gravitational force by buoyancy.In this case, the extension of the shoreline is limited by a steady-state value, which we identify as a function of the ratio of the densities of the ambient fluid and the current.When the submerged fluid exhibits an order-of-magnitude increase in viscosity, the behavior upon crossing the shoreline initially differs from the buoyancy case, with the formation of a high viscosity platform over which a subaerial current propagates.However, at sufficiently late times, the platform slumps, the material progresses as a subaqueous current, and the system approaches a steady shoreline position, as for the buoyancy-driven case (although with a greater extension than via buoyancy alone).In contrast, if, through hydroclastic fragmentation, the material delivered to the shoreline builds a static granular platform upon contact with the water, the shoreline can continue to extend indefinitely.We characterized the early and late time behavior of the current in this case, including the derivation of a similarity solution for the late time shape and rate of advance of an unconstrained delta, fed by a point source of material on the onshore slope.
We also explored buoyancy-driven inflation and platform building in analog laboratory experiments.Our experimental analog for the granular/viscous platform model involved an alginate solution which gels on reaction with calcium ions.In the absence of the reaction, the buoyancy-driven inflation of the current is predicted by the theory.On the other hand, when the bath contains calcium ions, the gelling reaction causes the current to progressively build a shelf over which the flow is diverted, thereby providing a qualitative analog for the platformbuilding models.Furthermore, the experiments show some interesting features including failure of the shelf near the original shoreline or at the tip of the shelf, and a subsequent propagation of a submerged current that is dynamically effected by a gel layer, but is nonetheless able to propagate as a coherent current.
In practice, of course, the behavior of lava upon entering water is highly complex, and our contribution here provides a first minimal mathematical model for the evolution of a lava delta.Further work could include modeling the role of seabed topography, accounting for temperature dependent rheology and the flow of lava through the porous granular platform, and determining the instability of the platform to gravitational collapse under the weight of the growing subaerial cap.which corresponds to the evolution equation for the granular platform model Equation 25, where x > 0. To resolve the sharp transition over x f < x < xb , we switch to a stretched coordinate ξ = ε 1 ( xb x) and write the leadingorder form of the evolution equation as where dot represents the dimensionless time derivative.Integrating, given ĥ = 0 at ξ = 0, we find But, because xb = x f + O(ε) and ĥ = x f at x = x f , the flux there is given by which must match with the flux at this location from x < x f .Explicitly, which is equivalent to the granular platform boundary condition, Equation 23, for the particular value of k = 1.Hence, the viscous platform model with ɛ → 0 and the granular platform model with k = 1 are equivalent.
For small but finite values of ɛ the models ultimately diverge from one another since the viscous platform model results in a steady-state delta while the granular platform grows indefinitely.We can determine the timescale on which the models diverge by examining the size of the neglected terms.For example, for the planar case, Equation A3 evaluated at ĥ = x = x f indicates where we previously neglected the second term in the brackets on account of being O(ɛ).However, for fixed ɛ and large times t, we have that x f ∼ ̅̅̅̅ 2 t √ (see Equation 29) from the late time behavior of the granular platform model.Thus the neglected term becomes O(1) and the equivalence breaks down when t ∼ (ε(1 R)) 2/ 3 .Figure A1 shows a comparison of depth profiles for the viscous platform model with ɛ = 2 × 10 3 and R = 0.7 and the granular platform model with k = 1, showing excellent agreement as anticipated by the analysis above.For the final snapshot, which is approximately 15 dimensionless time units after the current meets the water, the slight discrepancy between the models is shown.From Equation A6 with t = 15, we anticipate that xb x f ≈ 0.18, matching the order of magnitude of the discrepancy visible in the close-up shown in Figure A1.

Appendix B: Granular Platform: Derivation of Boundary Condition
Material from the surface lava flow enters the ocean and is assumed to transform rapidly into a granular media and to avalanche down the progressively growing submarine platform.This pile of granular material evolves quasistatically so that there is instantaneous dynamical balance between its submerged weight and the normal stresses, and between the streamwise gravitational forces and the basal resistance.Writing the dimensional stress tensor within the granular platform by σ ij , and the thickness of the granular pile above the seabed by h, vertical force balance requires that σ zz = ρ(1 R)g(h x tan α z). (B1) Then, on the assumptions that the normal stresses are isotropic and the upper surface of the submerged pile is stress-free, force balance in the horizontal directions yields where σ bx and σ by are the components of the basal shear stress aligned with the x and y axes, respectively.If the material is treated as Coulomb-like with coefficient of friction, tan θ r , we deduce that The elevation of the submerged platform must match the surface current at its boundary.The edge of the surface lava flow is located at (x, y) = (x s (t; s), y s (t; s)), where t denotes time and s the parameter that specifies the boundary.Thus h = x s tan α on x = x s and y = y s .(B4) Using Charpit's method (see, for example Ockendom et al., 2003), the solution to Equation B3 satisfying Equation B4 is given parametrically by h = x tan α 2r tan 2 θ r (B5) 1/2 , (B6) where the parameter r ≥ 0 and a prime denotes differentiation with respect to s.The characteristic curves, x(r; s) and y(r; s) given by Equations B6 and B7, and known as rays in this context, are straight lines that are perpendicular to the curve (x s , y s ) and outward orientated from the domain enclosed by the curve.It is convenient to eliminate the parameter r, in which case the solution is given by h = x tan α ((x x s ) 2 + y y s ) 2 ) 1/2 tanθ r .(B8) A particular aspect of this solution is that the maximum gradient of the submerged granular pile is the angle of repose, and this is attained along the direction specified by the rays Equation B6 and B7 (namely perpendicular to the evolving shoreline (x s , y s )).
At t = t 0 , the edge of the surface lava flow is at a location, S (see Figure 2d) with coordinates (x s (t 0 ; s), y s (t 0 ; s), 0).The seabed vertically below S is therefore at location, B = (x s (t 0 ; s), y s (t 0 ; s), x s (t 0 ; s) tan α).We define an angle, β, such that the unit normal to the curve (x s , y s ) is given by n = ( ) ≡ ( cos β, sin β).

(B9)
This means that the submarine platform extends to a location P with coordinates (x p , y p , x p tan α), where x p = x s (t 0 ; s) + 2r cos β tanθ r and y p = y s (t 0 ; s) 2r sin β tanθ r (B10) and since h = 0 at P, this implies that 2r tan 2 θ r = x p tan α.Finally we note that the angle PSB is π/2 θ r .The area of the triangle SBP, denoted by A, is then given by The submerged granular pile grows through the delivery of a volume flux of material at the edges of the surface current.The pile is assumed to adjust rapidly through granular avalanching to its quasi-static form derived above.Thus the volume flux delivered from the lava flow, Q, balances the rate of increase of A and the rate of advance of the platform, modulated by the solid fraction, ϕ.Then, writing the normal velocity u n = ẋs / cos β, where a dot denotes differentiation with respect to time, we find After non-dimensionalization, this gives the boundary condition

•
The dynamics of evolving lava-fed deltas are explored through shallowlayer mathematical models and analog laboratory experiments • The effects of reduced gravitational driving and formation of a viscous or granular subaqueous platform exhibit distinct late-time behaviors • Gelling of sodium alginate solution provides a novel laboratory analog for abrupt rheological changes which occur when lava enters water Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Example of lava delta formation on Kilauea volcano, Hawai'i.Photographed in May (left) and July (right) of 2017, showing the extension of the shoreline.Photograph used courtesy of the United States Geological Survey (DeSmither, 2017).

Figure 2 .
Figure 2. Sketches of the problem geometry for the water entry of a viscous gravity current.Panels (a)-(c) show the two-dimensional geometries (i.e., planar or Hele-Shaw), for the cases of (a) buoyancy-driven inflation, (b) the viscous-platform model, and (c) the granular-platform model.(d) Shows a schematic of an unconstrained, three-dimensional current (for which the three delta-building models can be defined analogously to (a)-(c)).The triangle indicated by points SBP is referred to in Appendix B when deriving the boundary condition at the nose of the subaerial flow.

Figure 3 .
Figure3.Profiles of the free surface of a planar, buoyancy-inflated current with R = 0.7, plotted at constant time intervals (a) as three-dimensional curves above the ( x, t) plane and (b) projected onto the x ẑ plane.In dimensionless coordinates, the underlying slope is located at ẑ = x.The blue line indicates the position of the shoreline ĥ = x) and the magenta line indicates the surface profile for the steady state Equation16.(c) Dimensionless position of the front (red) and shoreline (blue) as functions of dimensionless time, t.(d) Profiles of the free surface of a buoyancy-inflated current with R = 0.7 in the Hele-Shaw geometry, and the corresponding steady state prediction (magenta) Equation17, for comparison with (b).

Figure 4 .
Figure 4. Profiles of the free surface of a planar current under the viscous platform model with R = 0.7 and ɛ = 2 × 10 3 , plotted at constant time intervals (a) as three-dimensional curves above the ( x, t) -plane (at intervals of 30) and (b) projected onto the x-ẑ plane (at intervals of 10).In dimensionless coordinates, the underlying slope is located at ẑ = x.The blue line indicates the position of the shoreline ( ĥ = x) and the magenta line in (b) indicates the uniform subaqueous surface height for the steady state.The inset shows the surface profile (at constant time intervals of 4) prior to significant slumping of the viscous platform.(c) Dimensionless position of the front (red) and shoreline (blue) as functions of dimensionless time, t.The dashed blue line indicates the final, steady position of the shoreline, given by Equation22.(d) A close up of the gray rectangle in (c), which is prior to significant slumping of the viscous platform.

Figure 5 .
Figure 5. Profiles of the free surface of a planar viscous current propagating over an expanding granular platform with k = 0.64, plotted (a) as three-dimensional curves above the ( x, t) plane and (b) projected onto the x ẑ plane.The blue line indicates the position of the shoreline ĥ = x) and the magenta line indicates the surface profile for the steady-state submerged current from Figure 3. (c) Dimensionless position of the front (red) and shoreline (blue) for the solution shown in panels (a)-(b).In (d), the shoreline position, x f , is plotted logarithmically and over longer times; the red and blue dashed lines indicate the

Figure 7 .
Figure 7. Rheometric flow curves for the alginate and sugar solution.The symbols distinguish different samples and the up and down ramps of shear-rate controlled tests.The solid line shows the best-fit power law.The agreement between up and down ramps indicates that the fluid is not thixotropic.

Figure 6 .
Figure 6.Sketch of the experimental setup.

Figure 8 .
Figure 8. Sample images and extracted interfaces from three experiments of a gravity current of alginate and sugar solution entering an ambient bath.One experiment is shown for an ambient bath of DI water (a-c) and two are shown for a bath of CaCl 2 solution (d-f) and (g-i).The grid shows a scale of 1 cm.

Figure 9 Figure 9 .
Figure9shows the time series for the shoreline and front positions (panel a) and a sample free-surface profile (panel b) from the unreacting experiment in Figure6.Also shown is a solution of the model for buoyancy-driven inflation of a viscous current in a Hele-Shaw cell results from an additional unreacting experiment.To make this comparison, we estimate the volume flux per unit width, Q ≈ 0.088 cm 2 s 1 , by evaluating the rate of change of the area of red fluid visible in the video frame, and measure the incoming layer height and slope angle to be H ≈ 7 mm and tan α ≈ 0.2.Hence, L = H tan α ≈ 3.5 cm, T = HL/Q ≈ 27.8 s.(30)

Figure 10 .
Figure 10.Time series of dimensionless front ( xb , red) and shoreline ( x f , blue) positions of reacting alginate currents for four experiments.Dashed lines marked A indicate times at which the gel shelf failed close to the initial shoreline position, and dashed lines marked B indicate times at which the shelf failed at the leading edge of the shelf, resulting in a halting or retreat of the shoreline position.The dashed blue line shown in panel (b) shows the shoreline predicted by a numerical solution of the granular platform model detailed in Section 4.2 using the Hele-Shaw flux rule and k = 0.1.

Figure 11 .
Figure 11.Dimensionless current thickness, ĥ (color plots), of numerical solutions for three-dimensional viscous gravity currents entering water, for (a, b) buoyancy-driven inflation, and (c, d) viscous-platform building with ɛ = 10 3 .The entire computational domains are shown; R = 0.833.The red dashed line highlights the edge of the currents; the blue dashed line indicates the shorelines.(e) Cross sections through the line of symmetry, ŷ = 0, for the solutions in panels (b) (red) and (d) (blue).

Figure 12 .
Figure 12.(a) Similarity delta shape in similarity coordinates X = X s (Y) for three values of r = tan α/tan θ r (see legend).(b-d) Similarity delta shape for r = 0 (red dashed) against color plot of current thickness from full numerical simulation of the viscous platform model with ɛ = 10 3 .The numerical simulation snapshots are shown at dimensionless times of t = 75 (b), 150 (c) and 300 (d) from the initiation of the source at x = 4, while the similarity solutions are shown for ts = 44 (b), 100(c) and 220 (d), such that the total volume in the delta approximately matches the numerical simulation.

Figure 13 .
Figure 13.Sample images from experiments with unconstrained currents.Panels (a)-(c) show an unreacting experiment at 2 min intervals; panels (d)-(i) show an experiment with the alginate reaction, at 1 min intervals.The blue dashed lines locate the shorelines, and the red dashed lines in (a)-(c) locate the current edge, which is less clear in the unreacting experiment due to the absence of the sharp gel interface.In panels (j-l), side views of the evolving "delta" for the reacting current in panels (d-i) are shown at the indicated times.Panel (l) captures the appearance of a thin subaqueous current following a failure of the gel sack holding back the pool of unreacted alginate.All times are measured from when the current met the bath.

Figure A1 .
Figure A1.Snapshots of depth profiles for the viscous platform model with ɛ = 2 × 10 3 and R = 0.7 (black) and the granular platform model with k = 1 (red dashed).The snapshots shown are at dimensionless time intervals of 4, from the initiation of the unit flux at x = 5.The second panel shows a closeup of the front of the current at t = 20, where the discrepancy between the models is just beginning to be apparent.

Journal of Geophysical Research: Earth Surface
TAYLOR-WEST ET AL.
Together these indicate that