Some Lava Flows May Not Have Been as Thick as They Appear

Individual lava flows in flood basalt provinces are composed of sheet pāhoehoe lobes and the 10–100 m thick lobes are thought to form by inflation. Quantifying the emplacement history of these lobes can help infer the magnitude and temporal dynamics of prehistoric eruptions. Here we use a phase‐field model to describe solidification and remelting of sequentially emplaced lava lobes to explore additional processes that may lead to thick flows and lobes. We calibrate parameters using field measurements at Makaopuhi lava lake. We vary the lobe thicknesses and the time interval between eruptions to study the interplay between these factors and their impact on the thermal evolution of flows. Our analysis shows that if the time between emplacements is sufficiently short, remelting may merge sequentially emplaced lobes—making lava flows appear thicker than they actually were—which suggests that fused lobes could be another mechanism that creates apparently thick lava flows.

wide range of flow thicknesses. While our model could still be applicable to shallow flows undergoing predominantly unidirectional solidification and moving at low velocities, the model simulated in this paper cannot capture the inherently 3D, meandering nature of pāhoehoe flows, especially those present during the amalgamation of small ( 1 m), predominantly noninflated pāhoehoe lobes, for example, in lava piles where hummocky pāhoehoe or "compound" lavas are forming (e.g., Baloga et al., 2001;Hamilton et al., 2020). Thus, we focus most of our discussion on large CFB flow lobes where our model is best suited, since for such flows, the heat transfer in the nonvertical directions is dominated by that in the vertical (Hon et al., 1994;Wright & Okamura, 1977).
We use simplified magma solidification models to constrain how quickly two subsequent flow lobes must be emplaced to fully merge, thereby providing constraints on CFB eruption tempo. In Section 2, we describe a new phase-field model for lava lobe cooling, and then simulate the solidification of a single flow lobe and two sequentially emplaced flow lobes using our model in 1D. In Section 3, we use these results to outline three distinct regimes (fused, in parallel, in sequence) for interlobe cooling. Finally, in Section 4, we compare our results with observations to assess whether remelting can help explain thick CFB flows and analogous thick flows in other planetary settings. Our results are used to put lower bounds on how quickly CFB flow fields were emplaced in order to preserve multiple lobes within a single flow.

Model Equations
The phase-field framework is a mathematical approach to describe systems out of thermodynamic equilibrium (Anderson et al., 1998), first introduced in the context of solidification processes and phase transitions of pure or multicomponent materials (Boettinger et al., 2002;Cahn & Hilliard, 1958). The framework evolves the solidification front via a system of partial differential equations, avoiding the need for explicit tracking of the moving interface as traditionally done in the Stefan problem (Anderson et al., 1998). Here, we consider a simplified model of lava solidification where we track the binary solidification of lava through a phase variable, denoted ϕ (ϕ = 1 for the melt and ϕ = 0 for the solid phase), with corresponding temperature, T. The evolution of ϕ and T can be described by the following system of partial differential equations: where T m is the melting temperature of the lava, α is the thermal diffusivity, ω ϕ is the interfacial width coefficient, τ is the characteristic time of solidification (not the solidification time of a lobe), L is the latent heat of fusion for lava, and H is the energy barrier; see Table S1 in Supporting Information S1 for the values of these parameters used in this work, which are adopted from the typical thermal properties of basaltic melt (Audunsson & Levi, 1988;Cooper & Kohlstedt, 1982;Patrick et al., 2004;Peck et al., 1977;Worster et al., 1993;Wright & Marsh, 2016). ∇ is a partial differential operator defined in Text S1 in Supporting Information S1. Γ = Γ ( ) and Ψ = Ψ ( ) are auxiliary functions of the phase-field model (see Text S2 in Supporting Information S1). Because we are working with a binary phase-field model, we do not model the so-called "mush zone" that exists in actual lavas (e.g., Wright & Marsh, 2016). Consequently, we combine the solidus and liquidus temperatures as T m = 1,070° C, which is within the range of reasonable values reported in the literature for the liquidus of typical basaltic magmas (e.g., Cashman & Marsh, 1988;Wright & Marsh, 2016).
We impose convective and radiative boundary conditions at the surface while fixing the temperature at the bottom of the domain. Moreover, we integrate our phase-field equations over a sufficiently large domain such that the lower boundary does not influence the temperature and phase during solidification (see Text S2 in Supporting Information S1). Because the horizontal dimensions (kilometers) are much larger than the vertical scale (meters) for the flows of interest here, we perform our simulations in a 1D vertical dimension. Consequently, the conductive heat transfer will be primarily in the vertical direction. We provide additional details regarding the numerical scheme we used in Text S3 in Supporting Information S1.

Model Validation and Limitations
The phase-field modeling parameters τ and ω ϕ are derived in terms of measurable quantities in Text S1 in Supporting Information S1 using the approach in Kim and Kim (2005), and then calibrated based on field data collected from Makaopuhi lava lake (Wright et al., 1976;Wright & Marsh, 2016;Wright & Okamura, 1977). The calibration results show that the model agrees with the lava lake data for a range of parameters (see Figure S1 and Table S1 in Supporting Information S1); we ultimately choose ω ϕ = 3.22 × 10 −1 m and τ = 2.90 × 10 6 s in our simulations, both of which are well within these ranges.
As an additional test, we use the calibrated parameters from Makaopuhi lava lake to simulate measurements of inflating pāhoehoe lava flows on the Kīlauea volcano in Hawai'i, taken from Hon et al. (1994). The results (Figures S2 and S3 in Supporting Information S1) show decent agreement at depths deeper than ∼10 cm below the cooling surface, although there is noticeable disagreement near the surface. These validation efforts demonstrate that, while our model robustly captures macroscopic cooling of lava across multiple data sets, it lacks accuracy in describing the temperature evolution in the uppermost section of cooling lava (∼10 cm). One explanation is that the bubbles and vesicles that accumulate near the lava's surface (Audunsson & Levi, 1988;Cashman & Kauahikaua, 1997;Self et al., 1998) tend to decrease α near the surface (Keszthelyi, 1994). We also neglect the temperature dependence of α (Jaupart & Mareschal, 2010). While it is possible to include these effects in the model, it would also introduce additional parameters that are challenging to calibrate, and would also require resolving multiphase physics and chemistry at the sub-centimeter spatial scale and subsecond timescale. The fact that we have neglected to include such effects adds an increasingly non-negligible degree of uncertainty to our results as the lobe size decreases, the degree of which should be explored in future studies. Nevertheless, the model presented here, although simplified, still captures the first-order effects of latent heat and thermal diffusion that dominate lava cooling while allowing us to simulate cooling processes spanning from seconds to years. From our validation tests above, we see that any surface effects appear to be negligible for depths below ∼10 cm. This is especially true for larger lava lobes, where the influences of latent heat release and diffusion on the velocity of the solidification front and evolution of temperature profiles are especially greater in both spatial and temporal extent than those which are due to the surface and near-crust phenomena aforementioned (Patrick et al., 2004;Wright et al., 1976;Wright & Marsh, 2016;Wright & Okamura, 1977).

Setup for Lava Cooling Simulations
We use the model to perform two types of simulations. We first simulate solidification of a single lava lobe of thickness h to obtain the total time t h it takes to reach complete solidification. The results are used to design the second set of simulations, where we simulate sequential emplacement of two lava lobes of equal thickness h, separated by a time period of t emp . We consider 17 different lobe thickness h, from 0.1 to 20 m, to explore the behaviors of both thin pāhoehoe lobes (<1 m), as seen in recent Kīlauea eruptions (Lundgren et al., 2019;USGS, 2019), and thick lobes (≫1 m), as seen in Columbia River Basalt Group and other CFBs (Self et al., 2021). For the sequential emplacement simulations, we scale t emp relative to t h and explore nine different emplacement intervals for each thickness: Here, t emp is sampled along multiple orders of magnitude in order to capture a wide range of cooling times. The lower bound for the emplacement intervals is based on typical pāhoehoe-type flows (Anderson et al., 1999;Hon et al., 1994) while the upper bound is provided by examples from flood basalt provinces (Self et al., 1996;.

Results
We perform a total of 153 simulations of the sequential emplacement of two lava lobes and identify three distinct qualitative regimes of interlobe solidification. These regimes can be delineated based on the ratio between t emp and the conductive time scale, the latter of which is more precisely described by t h , but approximated here by h 2 /α to be more physically interpretable (see Figure S4 in Supporting Information S1 for how well the approximation t h ∼ h 2 /α holds). Below, we describe each regime in detail with examples in Figure 1 for the case of h = 10 m lava lobes.
In sequence (t emp > 0.06 h 2 /α): The first lava lobe completely solidifies before the second lobe is emplaced ( In parallel (0.01 h 2 /α < t emp < 0.06 h 2 /α): As indicated by the narrowing of both black contours in the top plot and the decreasing melt thickness in the lower plot with time, both lava lobes solidify for overlapping time, but the interface between them does not remelt (Figure 1, middle). Because the bottom lobe is hot, the collective cooling of both lobes is slower than for in sequence flows, as indicated by the decrease in slope in Figure 1 (bottom middle).
Fused (0 < t emp < 0.01 h 2 /α): After emplacement, the solidified portion of the lower lava lobe remelts completely, after which both lobes combine to form a single, larger lobe. For early times, there are four solid-melt interfaces that correspond to the simultaneous solidification of two independent lobes. However, the two interior interfaces eventually disappear, which marks the merging of the two lobes. The remelting event is also evident when we track the total melt thickness over time (Figure 1, right). After the arrival of the second lobe (indicated by the red dot), the total melt thickness increases slightly at some point, corresponding to the remelting that caused a reduction in the solid fraction. Despite a monotonic loss of entropy over time after the second flow arrives, remelting can still occur, since some sensible heat is converted into latent heat. In the other two regimes, the melt thickness never increases after the arrival of the second lobe.
We compile the results from all the simulations into a regime diagram in Figure 2, which shows the combined control of individual lobe thickness and emplacement intervals on the interlobe solidification during sequential emplacement. We map the three regions of interlobe solidification, separated by two boundaries extrapolated from our results: t emp = 0.01 h 2 /α and t emp = 0.06 h 2 /α. These regimes and the boundaries that define them are universal for both thin and thick lobes.
The bottom four panels in Figure 2 also illustrate examples of lava flows that appear to have been emplaced in parallel or in sequence, as suggested by their distinct interlobe boundaries. These examples are also marked in the regime diagrams, where the vertical position of the marker corresponds to the minimum emplacement interval predicted by our model (e.g., t emp = 0.01 h 2 /α). The hexagonal marker corresponds to ∼10 cm thin lobes seen in the Kupaianaha flow field  that are predicted to have been emplaced at least ∼4 min apart. The square marker corresponds to ∼0.5 m thin lobes seen in Elephanta Caves (Deccan Traps) (Patel et al., 2020), and are predicted to have been emplaced at least ∼2 hr apart. The circular marker corresponds to ∼8 m thick lobespart of a single flow and seen in Rajahmundry Traps (Fendley et al., 2020)-that are predicted to have been emplaced at least ∼20 days apart. The star-shaped marker corresponds to ∼20 m thick lobes seen in Columbia River Basalts (Self et al., 2021) that are predicted to have been emplaced at least ∼4 months apart.

Discussion
There is a body of the literature that commonly assumes that even the thickest ( 40 m) CFB flows were formed by flow inflation (e.g., Anderson et al., 1999;Rader et al., 2017;Self et al., 1996Self et al., , 1998, inspired by observations of Hawaiian lava flows (Hon et al., 1994). However, our analysis suggests that even thick (30-40 m total height) flows could have arisen by fusing lobes together if eruption intervals are shorter than a month or two. One practical challenge in testing our proposed mechanism is the ability to identify fused flow boundaries in the field, since fusing would remove structures corresponding to the crusts of the two lobes. However, some relics of the originally distinct flows may remain, such as compositional differences Vye-Brown et al., 2013) and possibly structures indicative of fused flow crusts, such as multiple differential cooling zones and vesicle-rich horizons (see Text S6 in Supporting Information S1 and Figure 3). Moreover, vesicle-rich horizons are commonly interpreted as remnants of inflation , and so the presence of these alone may not be sufficient to distinguish between lobe inflation and fusion, or at least with our current understanding of how the observable characteristics of said horizons reflect their formation. In the two scatter plots-for which the tick mark spacing scales quadratically along the vertical axis, that is, with 2 emp rather than t emp , so that the trend is roughly linear with h-we omit the top two lines of points corresponding to t emp = 2 3 t h , 2 4 t h for sake of visual clarity. The example 20 m scale lobes in CRBG are taken at Palouse Falls, Washington, USA (Sheth, 2017). See Figures S10 and S11 in Supporting Information S1 for high-resolution versions of the last two examples.  a, 119°54′52″W, 46°47′50″N). Coll m -Columnar, Ent-Entablature, and Col-Colonnade-the jointing pattern is divided into three categories: Col, Ent, and Coll m Ent (which is intermediate between the two end-members). The right panels show the vesicle porosity and geochemical variations in the DC-16 borehole. Panels (c, d) show stratigraphic sections with geochemical and textural variations in the Cohassett flow in the RRl-6 Core and DC-6 cores respectively; data from . We also show the assigned compositional types to parts of the Cohassett flow by .

Potential Example of a Fused CFB Flow
One potential example of a fused CFB flow is the ∼70 m thick Cohassett Flow from the Columbia River Flood Basalts. The flow is part of the Grande Ronde Basalt Group and is a member of the Sentinel Bluffs Member lava flows in Pascoe Basin (e.g., McMillan et al., 1989;, see Figure 3a for a map of outcrops and drill core data). As shown by the annotated picture in Figure 3b, the Cohassett has a multitiered structure with alternating entablatures and colonnades (see Text S6 in Supporting Information S1 for a description), as well as a 6.5 m thick internal vesicular zone (IVZ; ∼20 m from the flow top, Figures 3b-3d) with many ∼1 cm diameter vesicles (Mc-Millan et al., 1989;Tomkeieff, 1940). To first order, the Cohassett flow in the outcrops (Figure 3) appears to be a single thick sheet lobe. The Cohassett flow also exhibits one of the most striking geochemical variations amongst the Grande Ronde flows. The flow has an approximate vertical bilateral symmetry geochemically centered just under the IVZ, as seen from data across sections more than 50 km apart (Figure 3). Using characteristic patterns in TiO 2 , P 2 O 5 (and other major and trace elements),  defined four distinct compositional types within the flow: California Creek, Airway Heights, Stember Creek, and Spokane Falls. Typically, these compositional types are separated by a vesicular horizon. For example, a horizon ∼13-15 m from flow top separates massive basalt of the California Creek composition from the Airway Heights composition. Similarly, the Airway Heights and Stember Creek transition is characterized physically by a series of large vugs. The IVZ acts as the contact between the Spokane Falls and the Stember Creek compositional types (Figures 3b-3d). Finally, a vesicular horizon ∼40 m from flow top defines the transition from the Spokane Falls back to the Stember Creek compositional types. Interestingly, the subsequent compositional type changes from Stember Creek to California Creek/Airway Heights lack clear vesicular horizons (Figure 3).
Corresponding spatially with these geochemical changes, the Cohassett flow also exhibits systematic changes in plagioclase abundance and fine-grained fraction (groundmass, Figure 3c based on data from Reidel, 2006). In particular, the flow part comprising the IVZ and the Spokane Falls composition member has a fine fraction significantly more indicative of a flow top rather than the flow interior. Thus, this flow interior was potentially emplaced rapidly and cooled faster than a continuously inflating flow lobe interior (McMillan et al., 1989;Philpotts & Philpotts, 2005). The IVZ-entablature-colonnade sequence in the Spokane Falls lava further supports the conclusion that the cooling rates in this part of the flow were more akin to a flow top (DeGraff et al., 1989;Forbes et al., 2014). Even on an overall flow scale, the textural data for the Cohassett flow are inconsistent with the slow cooling expected for a ∼70 m flow; the plagioclase crystal size does not significantly change throughout the flow, unlike the case for a slowly cooling ponded lava lake (Cashman & Marsh, 1988;Philpotts & Philpotts, 2005).
Previously,  proposed that the Cohassett flow formed via the combination of different sheet flows (for each compositional type), each sourced from a different magma reservoir and eruptive vent. These individual flows sequentially intruded into the Cohassett flow as flow lobes and inflated it to its final height. However, the  model does not explain the abrupt shift to distinct compositional types along with sharp vesicle horizons (Figures 3b and 3b1-b2) without any signs of magma mixing or shear instabilities, despite intrusion and transport within the Cohassett flow for 10s of km. Alternatively, Thor Thordarson (personal communication; see also Vye-Brown et al., 2013) proposed that the Cohassett flow was formed by semicontinuous inflation with changing magma compositions in the magmatic system feeding the eruption. As evidenced by observations from some modern long-lived basaltic eruptions, for example, Pu 'u 'Ō 'ō eruption at Kīlauea, Hawai'i from 1983 to 2018 (Garcia et al., 2021), these changes can be relatively abrupt and could correspond with the presence of a vesicle horizon. Philpotts and Philpotts (2005) proposed that crystal-mush compaction in an inflated sheet lobe can also partially explain the observed geochemical variation and bubble segregation.
Here, we put forward a third alternative, building upon the original idea proposed by : We posit that the Cohassett flow is an example of a fused flow with multiple flow lobes. Suppose the Cohassett was close to the boundary between the fused and in parallel flow types (Figure 2). In that case, the presence of separating vesicle horizons as well as high fine-grained size fraction, especially for the Spokane Falls type, can be explained. Within this scenario, each constituent ∼10-20 m lobe would have to be emplaced within a few months of the previous lobe. However, more detailed modeling work specifically focused on the Cohassett as well as textural analysis, for example, stratigraphic crystal size distributions to estimate cooling rates (Cashman & Marsh, 1988;Giuliani et al., 2020), would be needed to ascertain which of the proposed models is correct and if Cohassett is indeed a fused flow.
It is similarly challenging to distinguish between in parallel and in sequence flows based on field volcanological observations alone without detailed textural analysis. One potential distinguishing feature may be the 2D shape of the bottom flow lobe in an in parallel flow since it will be viscoelastically deformed by the load from the overlying flow lobe (Abbott & Richards, 2020). One consequence of this would be the formation of squeeze-up structures at flow lobe edges seen in some CFB flow edges, for example, for the Western Ghats and the Rajahmundry Trap flows in the Deccan CFB Fendley et al., 2020).

Relevance of Fused Flows for Planetary Geology
Our results also have implications for inferring eruption conditions on other planetary bodies (Venus, Mars, Mercury, the Moon) where we can only observe the final lava flow thickness from remote sensing observations. Below, we briefly summarize observations of thick lava flows on each of these bodies: 1. Venus. With a surface area (∼60,000 km 2 ) comparable to many CFB flows (Moore et al., 1992;Wroblewski et al., 2019), the most striking example of a thick lava flow (∼100-200 m) on Venus is the Ovda Flactus flow. Although initially classified as a potentially high silica or rhyolitic flow, morphological analysis by Wroblewski et al. (2019) shows that Ovda Flactus has an emplacement rheology most consistent with basalt. Several other basaltic flows mapped across Venus have thicknesses ranging between 30 and 100 m (Guest et al., 1992;Lancaster et al., 1995;MacLellan et al., 2021;Moore et al., 1992;Zimbelman, 2003). We envision that our proposed process of fused flows can help explain these large flow thicknesses, especially since cooling rate during solidification on Venus will be smaller (∼30-40% compared to Earth) due to higher surface temperatures (Snyder, 2002). Correspondingly, individual flow lobes can be separated by longer times and still fuse together 2. Mars. There is a wide range of estimated flow thickness for Martian lava flows with values ranging from a few meters to ∼100-125 m (Hiesinger et al., 2007;Mouginis-Mark & Rowland, 2008;Mouginis-Mark & Yoshioka, 1998;Peters, 2020;Zimbelman, 1998). In particular, both the Tharsis volcanic province and the Elysium Planitia region on Mars have a number of flows with typical thickness greater than 40 m (Peters, 2020). These flow thicknesses are challenging to explain with inflation, but may be explained by lobe fusion. Since the mean Martian surface temperature (−60 °C) is colder than Earth's, the lava solidification time is about 5% shorter than on Earth (based on results from a model with Martian surface temperature). Since Martian surface gravity is ∼38% of Earth's, flow lobes may inflate to greater thickness before overcoming basalt yield strength to form new breakouts 3. Mercury. Du et al. (2020) used observations of partially and completely buried impact craters on Mercury to estimate lava flow thicknesses and found values between 23 and 536 m with a median of 228 m (with potentially even thicker flows), consistent with some other estimates, for example, 180 m (Wilson & Head, 2008). Even if we account for the difference in surface gravity (38% of Earth's), the median thickness of 228 m translates to 86 m of Earth's equivalent (in terms of flow dynamics; the cooling times are independent of gravity) 4. Moon. Lunar mare basalts have a range of flow thicknesses, from thin ( 1 m) to thick ( 100 m). These estimates come from a combination of in situ observations by Apollo Astronauts, remote sensing, and lunar penetrating radar on rovers (Chen et al., 2018;Gifford & El-Baz, 1981;Hiesinger et al., 2011;Rumpf et al., 2020;Spudis & Guest, 1988). Since lunar gravity is only about 16% of the terrestrial gravity, lobe inflation could form much thicker flows (a 15 m flow lobe on Earth will be equivalent to a 91 m lobe on Moon with respect to lateral pressure gradients). Thus, except for the thickest flows ( 100 m), flow fusion may not be required to explain the observed lunar flow thicknesses In aggregate, there are a number of very thick planetary basaltic flows. We posit that these could be potentially emplaced by the same process as we are proposing for thick terrestrial CFB flows (fused flow lobes). However, more data and careful analysis is necessary to rule out the possibility that these were primarily formed via flow lobe inflation.

Implications for Eruption Rates
In combination with viscous flow models for channelized lava flows (e.g., Jeffrey's equation), measurements of lava flow thickness have been used to estimate the mean flow velocity and viscosity (e.g., Baloga et al., 2001Baloga et al., , 2003 flow velocities are used in combination with constraints on total flow volume to calculate an eruption duration. However, these calculations are predicated on the assumption that the final flow thickness is representative of the molten channelized flow thickness at the time of eruption. Since the velocity depends strongly on the lava flow thickness (velocity ∼ thickness 2 for a Newtonian fluid), an incorrect estimate could strongly impact the estimates of eruption rate (∼flow rate × thickness) . Suppose that the observed lava flow results from the merging of two (or more) equally thick flow lobes. Then, the models will overestimate instantaneous effusion rates by about a factor of 2 3 = 8 (or larger for three or more fused lobes) if the final flow thickness is used as the characteristic open channel lava flow thickness, that is, with an error that grows ∼ cubically. This issue is further exacerbated by how, during the emplacement of pāhoehoe flows by inflation, the "active" molten part of the flow is smaller than the thickness of each lobe.
In the other end-member, wherein the total observed flow thickness is assumed to be only a consequence of lava flow inflation, most models require a relatively continuous, long-lived effusion rate that gradually thickens the flow lobe (Hamilton et al., 2020;Hon et al., 1994;Kauahikaua et al., 1998;Self et al., 1998). However, if the CFB flow is a fused flow, each constituent flow lobe must inflate to a smaller thickness. Consequently, since the two lobes can be sequentially emplaced days to months apart, the total eruptive duration can be smaller and/or allow for more effusion rate variations but will still form a single fused flow in the end. In conclusion, the possibility that different flow lobes merged into a single flow has important implications for inferences of effusion rate, especially its steadiness.

Conclusion
We provide a theoretical lower bound on emplacement interval that distinguishes a fused flow from nonmerged flows. For instance, a distinct boundary between two lobes of 10 cm each suggests that they are emplaced at least 4 min apart (t emp > 0.01 h 2 /α ≈ 4 min). The same calculation for two 20 m thick lobes suggests that the emplacement interval is at least 4 months if a distinct boundary exists between the two lobes. Furthermore, while it is often assumed that the 10-100 m thick lobes found in flood basalt provinces are primarily formed from inflation, our results suggest that these large lobes could also have been formed by smaller lobes emplaced in quick succession. Relatedly, some volcanological studies could be overestimating eruption rates and flow velocities by assuming that solidified flows belong to one flow rather than a series of smaller flows that merged with little to no trace of their original separation. While more work is necessary (especially empirically) to distinguish conclusively between a large, solidified flow that formed primarily via multiple-lobe emplacement versus the usually assumed mechanism of inflation, we still propose that the emplacement and subsequent fusion of multiple lobes is a plausible, additional process for forming CFB flows.
Noting some limitations, we also demonstrate the effectiveness of using phase-field models in simulating observed lava solidification over a range of timescales faithfully (Hon et al., 1994;Wright et al., 1976;Wright & Marsh, 2016;Wright & Okamura, 1977), with some local deviation ∼10 cm near the top surface. The phase-field model can be generalized to arbitrary domains in higher dimensions and account for additional complexities in thermal diffusivity, flow, and nonlinear rheology.

Data Availability Statement
All relevant simulation data, movies, figures, and codes can be found at https://zenodo.org/badge/latestdoi/357729300. In particular, all data files are contained in finaldata.zip, all codes can be found in the folder finalcodes, all relevant figures from this paper can be found in the folder finalfigures, and all movies (plus some extra movies) can be found in the folder movies.