Lunar Mare Lava Flow Dynamics and Emplacement: Predictions of Non-Newtonian Flow Dynamics, Syn- and Post-emplacement Cooling and Volatile Release Patterns, and Vertical and Lateral Flow Structure Development

We apply basic principles of magma ascent from deep source regions and its eruption into a low-gravity vacuum environment to develop a theoretical treatment of the fluid dynamics and thermodynamics of mare basalt lava flow emplacement and evolution on the Moon. The vacuum conditions influenced the release of volatiles in magma passing through lava fountains, thus controlling the syn- and post-emplacement vesicularity of the resulting deposits. To explain observed lengths and volumes of Mare Imbrium–type flows, high (106–105 m3 s−1) initial magma eruption rates were needed. Combined with low lunar magma viscosity, these caused flows to be initially turbulent. Resulting high radiative heat loss and consequent high crystallization rates caused rapid non-Newtonian rheological evolution and suppression of turbulence at tens of kilometers from vents. Slower cooling rates in the subsequent laminar parts of flows imply distinctive crystal growth rate histories. In a four-phase sequence, (i) initial transient dike-tip gas release followed by (ii) Hawaiian fire fountain activity with efficient volatile loss (iii) transitioned to (iv) Strombolian explosions in a lava lake. Late-stage lava now able to retain volatiles intruded and inflated existing flow deposits after flow front advance ceased. Volatiles forced out of solution by second boiling as lava cooled caused additional inflation. Low gravity and lack of atmospheric pressure commonly produced very vesicular lava. Escape of such lava through cracks in flow crusts is a possible source of ring-moat dome structures; collapse of such lava may explain irregular mare patches.


Introduction
Lava flows on Earth emerge in subaerial, submarine, and subglacial environments (Figures 1(a), (b); e.g., Lockwood & Hazlett 2010;Harris & Rowland 2015;Houghton et al. 2015) and have been used as a baseline for analyzing and interpreting lava flows on other planetary bodies (Wilson & Head 1983), including the Moon (Figure 1(c); Wilson & Head 1981).Basaltic eruptions in the subaerial terrestrial environment are typically characterized by a series of phases ("curtain of fire" fissures, vent centralization, Hawaiian activity, effusive lava flows, Strombolian activity) depending on their total volume, volume flux, and volatile content (e.g., Harris & Rowland 2015).Key controls on terrestrial subaerial eruptions of magma with a given composition and volatile content are Earth's gravity and resulting overburden lithostatic pressure, and the presence of a 1 bar atmosphere (Figures 1(a), (b), (d)).Both factors influence patterns of volatile release and pyroclast dispersal (Figure 1(d)).Also important is the geologic environment: Shallow magma reservoirs propagate dikes laterally as in Hawai'i (e.g., Ryan et al. 1981) and mantle diapir source regions feed large igneous province flood basalt eruptions (e.g., Coffin & Eldholm 1992, 1994;Ernst & Buchan 1997, 2001;Ernst 2014).These basic factors dictate many of the characteristics of terrestrial subaerial lava flow surface vents and flow morphology (such as pāhoehoe and 'a'ā surface texture), flow characteristics in map view (levees, channels, etc.; Figure 1(e)), flow cross-sectional shapes, and detailed internal structure (e.g., presence, size, shape, and abundance of vesicles, colonnades, and entabulature; Figure 1(f)).
In contrast to the Earth (Figures 1(a), (b)), the Moon (Figure 1(c)) is characterized by the lack of an atmosphere, low acceleration due to gravity, low density of an anorthositic crust overlying the lunar mantle, and the lack of divergent and convergent plate boundaries (Hiesinger & Head 2006).These environmental factors greatly influenced the physics of the ancient explosive and effusive eruptions of the basalts that dominated lunar volcanism (Wilson & Head 1981, 2017a;Head & Wilson 1992, 2017;Morgan et al. 2021).Determining the morphology and structure of lava flows commonly more than ∼3 billion years old is made difficult by the presence of an overlying several meters of impact-produced regolith (McKay et al. 1991), which can develop to depths greater than the thickness of the original flow (Wilcox et al. 2005;Bart et al. 2011;Prieur et al. 2018).Remote-sensing data have successfully identified and mapped a series of basaltic lava units on the basis of their mineralogy interpreted from spectral reflectance (e.g., Pieters & McCord 1976;Hiesinger et al. 2002Hiesinger et al. , 2011)).However, only in a few cases have units been documented that plausibly represent discrete individual lava flows whose boundaries can be identified (Schaber 1973) and mapped in detail (Bugiolacchi & Guest 2008;Morota et al. 2011;Zhang et al. 2016aZhang et al. , 2016b;;Morgan et al. 2016;Chen et al. 2018).Exposed sections through mare lava flows outcrop only on steep slopes, where accumulating regolith is shed downslope, and bedrock is exposed in only a few places, such as in the walls of sinuous and linear rilles (Howard et al. 1972;Basilevsky et al. 1977;Florenskii et al. 1978), collapse pits and skylights (e.g., Haruyama et al. 2009;Robinson et al. 2012;Wagner & Robinson 2014), and impact crater rim crests and  Aubele et al. (1988), Self et al. (1996Self et al. ( , 1997)), and Thordarson & Self (1998).
crater walls.Thus, there is little outcrop information on the details of lava flow emplacement in terms of gas-release patterns, relation to explosive activity, flow surface morphology, crosssectional structure, and how these vary in time and space.In this contribution, we formulate a theoretical model of the eruption of magma on the Moon, emphasizing the nature and evolution of gas-release patterns, the fluid mechanics of lava flow emplacement, and the structure and morphology of flows during and immediately following their eruption.Predictions derived from this model can be utilized to suggest and interpret orbital and in situ observations and measurements that can test and refine the model, and also to assess the development of regolith on the lava flows following their cooling and solidification (e.g., Kerber et al. 2018).

Lunar Basaltic Eruptions
Despite the lack of detailed information on the internal structure of lunar lava flows, decades of orbital imaging and remote sensing have provided a firm basis for theoretical analysis and predictions concerning the generation, ascent, and eruption of basaltic magma on the Moon (Basaltic Volcanism Study Project 1981; Wilson & Head 1981, 2017a, 2018;Head & Wilson 1992, 2017;Jolliff et al. 2006;Morgan et al. 2021).Observations of the products of lunar volcanic eruptions (see summary in Head & Wilson 2017) provide evidence supporting these theoretical expectations.
Specifically, volcanic eruptions on the Moon were controlled by the interplay between two main factors that vary with time during an eruption: (i) the release of volatiles, and (ii) the dense rock equivalent volume release rate (flux) of magma from the vent (Wilson & Head 2018).Modern technology applied to Apollo samples has revealed that lunar magmas transferred a range of volatiles from the mantle (e.g., Saal et al. 2008;Hauri et al. 2011).The picritic magmas sampled as pyroclastic glass beads by the Apollo missions released up to ∼1000 ppm of carbon monoxide and up to ∼1000 ppm water, with perhaps ∼400 ppm of sulfur compounds (Rutherford et al. 2017).Other types of basalt may have contained smaller amounts of volatiles, especially water; values of ∼100 ppm have been estimated for KREEP basalts (Robinson et al. 2016), but Milliken & Li (2017) report spectroscopic detection of up to 400 ppm residual water in pyroclastic deposits, implying that the pre-eruption value must have been higher.Carbon monoxide was largely released at pressures corresponding to depths greater than 50 km (Figure 1(d)) whereas the remaining volatiles were predominantly released within 500 m of the surface (Rutherford et al. 2017).Due to the lack of an atmosphere, even very small amounts of volatiles have the potential to cause explosive eruptions in all lunar mafic magmas as they approach the surface (Figures 1(d), 2; Wilson & Head 1981), and we infer that almost all lunar lava flows were produced by eruptions accompanied by some form of explosive activity (Wilson & Head 2018; see Figure 2).
Magma reached the lunar surface in large volumes through a crust in which it was commonly very negatively buoyant (Figure 1(c)) with minimal physical or chemical interaction (Finnila et al. 1994).One mechanism allowing this involves the formation of very large dikes that grow slowly from partial melt zones at depths of several hundred kilometers in the mantle.The dikes disconnect from their sources, migrate upward very rapidly due to their positive buoyancy in the mantle, and stall centered on the crust-mantle boundary as they reach neutral buoyancy (Wilson & Head 2017a).If such dikes have a great enough vertical extent, their tops penetrate through the crust, producing a fissure vent at the surface.Alternatively, a laccolith-like magma body may be formed at a neutral density level by the accumulation of melt from multiple dikes rising from deep sources.A similar body could form by melt segregation into the upper part of a diapiric body stalled at a rheological trap at the base of the lithosphere (Head & Wilson 2017).Arrival of an unusually voluminous dike overpressurizes the laccolith and drives a dike to the surface.
Whatever its formation mechanism, the dike feeds an eruption after quickly losing the gas that has accumulated in the lowpressure region in the dike tip (Phase 1 in Figures 2(a), (b)).Based on considerations of dike size, geometry, and interaction with the lithosphere, Wilson & Head (2018) predicted that the initial volume eruption rate is large, ∼10 6 m 3 s −1 , while the lower part of the dike is still rising toward equilibrium (Phase 2 in Figures 2(a), (c)), decreases to ∼10 5 m 3 s −1 over the course of a few days as the dike approaches equilibrium, and then decreases further, to ∼10 4 m 3 s −1 , over the course of about a day as buoyancy equilibrium is reached (Phase 3 in Figures 2(a), (d)).Subsequently, lithospheric forces, compressive during the latter half of lunar geologic history (Solomon & Head 1980), cause the dike to become narrower (Phase 4 in Figures 2(a), (e)).The volume of magma erupted in both the early, pre-equilibrium phase and the late, post-equilibrium phase depends on the vertical extent of the dike relative to the thickness of the crust; the total magma volume is predicted to range from a few hundred cubic kilometers to ∼1000 km 3 (Wilson & Head 2017a).
Estimating the volumes of lunar lava flows requires information on their thicknesses and surface areas.Areas are readily obtained from images when flow boundaries can be identified, and thicknesses have been estimated by various methods.Radar sounding from landed spacecraft (Yuan et al. 2017;Chen et al. 2018;Lai et al. 2020) is able to penetrate though multiple subsurface layers, but since it is most sensitive to major discontinuities in physical properties and becomes less sensitive with increasing depth, it may combine flows separated by very thin regolith layers.Similarly, partial burial and degradation of impact craters (Du et al. 2019) can provide estimates of the total thickness of lavas emplaced after the formation of a crater, but again this may be the result of more than one flow embaying the crater.Measurements of the depths of craters that have penetrated through a surface layer into a spectrally distinct buried layer also allow estimates of the upper layer thickness (Hiesinger et al. 2002).Measurements using shadow lengths of the 30-35 m heights of some flow fronts in Mare Imbrium were made by Schaber (1973), but examples of such well-preserved flow fronts are not common.Measurements of the vertical extents of exposures in pit crater (Robinson et al. 2012) and impact crater (Rumpf et al. 2020) walls offer the chance to measure multiple flow sequences.The thicknesses of eight layers identified in the walls of the lava pit in Mare Tranquillitatis were interpreted to imply flow thicknesses from 3 to 14 m by Robinson et al. (2012) and the average thicknesses of layers visible in the walls of the craters Kepler, Dawes, Bessel, and Lichtenberg measured by Rumpf et al. (2020) were 6.6, 5.6, 9.7, and 7.7 m, respectively.Recent predictions of the expected wide range of internal vesicularity structures of lunar lava flows (Wilson et al. 2019) and the ways in which these respond to regolith generation by impacts (Head & Wilson 2020) suggest that large contrasts in material properties can exist between the interiors, basal layers, and tops of flows.Finite resolution in current images means that fragmental tops, vesicular interiors, and massive basal parts of flows may be readily misidentified as separate units, leading to significant underestimates of flow thicknesses.The above observational thickness estimates should probably be doubled, and thus imply typical flow unit thicknesses of ∼10 to 20 m.Thicknesses of this order, coupled with surface areas of individual large flows like those in Mare Imbrium, with lengths >300 km and widths averaging ∼40 km (Schaber 1973), yield volumes of ∼200 km 3 , consistent with the theoretical expectations of Wilson & Head (2017a) discussed in the previous paragraph.
Figure 2 summarizes the four phases of behavior predicted for typical lunar eruptions (Wilson & Head 2018).In this contribution these predictions are used to develop a systematic model of magma eruption, volatile release, initial lava flow emplacement, and subsequent development.We then compare the model results with the observations of flow morphologies described above to refine eruption rates and durations.

General Issues
The ranges of dense rock equivalent magma volume eruption rate, F, and magma volatile content, n, exert a critical role in determining the nature of the lava flows created by lunar eruptions (Figure 2).We first consider the influence of F on eruption style.The volume flux is the product of the dike length and width and the magma rise speed.A large dike width minimizes friction and produces a large rise speed.It follows that if the volume flux is large, it is likely that the dike width and the rise speed of the magma in the dike are both large.The magma rise speed controls the transit time of the magma from depth.As magma rises it releases volatiles as gas bubbles that nucleate with diameters of order 10 μm (Sparks 1978;Le Gall & Pichavant 2016), and these are much less dense than the surrounding liquid and rise through it (Figures 2(b), (c)).The initially small bubbles rise very slowly through the liquid magma, but as earlyforming bubbles rise to shallower depth the pressure reduction causes them to expand and accelerate.In eruptions on Earth, small new bubbles are constantly forming, and the potential exists for larger, faster bubbles to overtake and coalesce with smaller bubbles.The resulting increase in size causes further acceleration.If the magma rise speed is small enough, so that a large amount of time is available for this process to occur, a runaway situation can be reached in which giant bubbles, called slugs when they fill most of the width of the dike, rise to burst intermittently through the magma at the surface causing Strombolian explosions (Head & Wilson 2017, their Figure 10;Figure 2(d)).Recent analysis of lunar eruptions (Rutherford et al. 2017) suggests that CO production occurred over a restricted range of depths somewhat greater than ∼50 km (Figures 1(d), 2).This maximized the time available for relative bubble motion.Wilson & Head (1981, their Figure 13) estimated that for magma volatile contents of 250, 750, and 2000 ppm CO, the mean magma rise speeds that must be exceeded to ensure steady Hawaiian (Figure 2(c)) rather than intermittent Strombolian activity were ∼0.2, 0.3, and 0.5 m s −1 .However, their calculation was based on an early estimate, 3.5 km (Sato 1976(Sato , 1977)), of the depth at which CO production began.Rutherford et al. (2017) revised the depth to ∼50 km and suggested a maximum CO content of ∼1400 ppm.Since it is the magma transit time that determines the efficiency of gas bubble coalescence, the new depth estimate can be used to update the critical magma rise speeds: For volatile contents of 250, 750, and 1400 ppm CO they are 2.9, 4.3, and 5.0 m s −1 , respectively.These values are consistent with the estimates by Wilson & Head (2018) of the magma rise speeds typifying the transition between Phases 2 and 3 of typical lunar eruptions, as shown in Figure 2(a).
In contrast, if the magma rise speed is large enough, coalescence is minimized and a nearly uniform spatial distribution of bubbles, expanding from their initial nucleation diameters of ∼10 μm, causes fragmentation of the magma into a nearsteady stream of released gas and pyroclasts in a Hawaiian lava fountain eruption (Head & Wilson 2017, their Figure 9(a)).In the lunar case, the extreme expansion of the released gas into the lunar vacuum (Figures 1(d); 2(c)) disrupted the liquid magma into pyroclastic droplets with diameters ranging from a few tens to several hundred micrometers, as sampled by the Apollo missions (Scott et al. 1972;Schmitt & Cernan 1973;Arndt et al. 1984;Arndt & von Engelhardt 1987) and predicted theoretically (Wilson & Head 1981;Morgan et al. 2021).For the most volatile-rich magmas, the expansion of the released gas would have accelerated monodisperse (similar particle size) pyroclasts to speeds of ∼130 m s −1 , leading to ballistic ranges of ∼10 km.The range would have increased to ∼30 km for the fine size fraction if the particle size distribution was extremely bimodal (two well-separated peaks in the size distribution; Wilson & Head 2017a).Morgan et al. (2021) have recently shown that the distribution was probably bimodal.However, they also found that the decoupling process between pyroclasts and gas, as the gas expands to pressures so low that the normal gas laws do not apply, is more complex than previously assumed.This leads to predicted maximum ranges, R, of ∼20 km in steady lunar Hawaiian eruptions (Figure 2(c)).
The mean magma rise speed in the dike, U d , can be linked to the erupted volume flux, F, as follows (Figure 2(a)).At depths greater than ∼10 km the volume occupied by CO released from the magma is less than 10% even for n as large as 2000 ppm, and so the bulk magma density is essentially equal to the vesicle-free liquid density.At these and greater depths the relationship between U d and F is therefore given by the continuity requirement, where W and L are the mean width and the horizontal length of the dike.We do not expect the horizontal length of the dike to change greatly between the surface and depths of order 10 km, as eruption simulations for lunar conditions (Wilson & Head 1981, 2017a) show that the decrease in bulk magma density caused by gas expansion as the pressure decreases is largely compensated by the increase in mean magma rise speed due to that same gas expansion.We therefore also identify L as the length of the active fissure vent as the dike reaches the surface.Few examples of fissure vents have been identified in the lunar maria.Head & Wilson (2017) attribute this to the likelihood of fluid lava draining into and drowning the vent at the end of an eruption.Nevertheless, they give five examples (their Figure 24) of features such as lines of small cones or ridges that may mark the location of fissure vents, and these have outcrop lengths ranging from 9.1 to 13.0 km, with an average of 11.1 km.Also, the proximal sections of well-preserved lunar lava flows (Bugiolacchi & Guest 2008;Morota et al. 2011;Morgan et al. 2016;Zhang et al. 2016aZhang et al. , 2016b;;Chen et al. 2018) have widths of 10-30 km, consistent with some lateral spreading of lava as it travels downslope away from a fissure that is initially ∼10 km long, and we adopt this value for quantitative illustrations of a number of processes in what follows.Wilson & Head (2017a, their Table 2) explored the geometries of lunar dikes feeding surface eruptions and found that the ratio of total dike length to mean dike width ranged from 1050 to 1670 with an average of 1315.We adopt this average ratio of L/W together with a range of values of L from 5 to 15 km in Equation (1) to show how U d varies with L for a series of values of F in Figure 3.The shaded area marks the zone between magma rise speeds of ∼0.5 m s −1 (for CO-poor magmas) and ∼5 m s −1 (for CO-rich magmas), where we expect the transition from Hawaiian activity (at higher speeds) to Strombolian activity (at lower speeds) to occur (Figure 2(a)).For our nominal fissure length of 10 km, all eruptions with F greater than ∼3 × 10 5 m 3 s −1 are expected to produce steady Hawaiian fire fountains, and all with F less than ∼3 × 10 4 m 3 s −1 are expected to exhibit intermittent Strombolian activity, with increasingly pulsating fire fountains forming as the eruption rate decreases through this range.We now consider how each of the four phases of an eruption (Figure 2) influences the lava-forming processes.

Phase 1
This phase lasts for only a few minutes after the dike arrives at the surface (Figures 2(a), (b)).Any free gas at the upper dike tip (predominantly CO) is vented first.It expands from a pressure of ∼0.5 MPa (Wilson & Head 2017a) and reaches a speed in excess of 1 km s −1 , so that any small regolith particles that might be entrained in the gas are deposited over a large fraction of the entire lunar surface, forming an entirely undetectable deposit.This gas release is followed immediately by the explosive eruption of the extremely vesicular magmatic foam that has accumulated just below the free gas.This material represents only about 0.25% of the total mass of magma in the dike (Wilson & Head 2017a) and is dispersed over a wide area (∼10-30 km radius), forming a very thin deposit of pyroclasts (a few centimeters) that is likely to be buried by subsequent eruptive activity or mechanically mixed into the regolith as the latter gets thicker with time (e.g., Head & Wilson 2020).

Phase 2
This is the earliest significant phase of the eruption, characterized by Hawaiian activity (Figures 2(a), (c)).After their acceleration in, and just above, the vent (Morgan et al. 2021), the submillimeter-sized droplets into which magma is fragmented decouple from the expanding gas and follow ballistic trajectories forming a lava fountain.Pyroclasts fall to the surface to accumulate into a lava lake around the vent if they remain hot or to form a pyroclastic cone if they have cooled sufficiently.In practice, the typical sizes of lunar pyroclasts (∼10 μm to a few millimeters: Heiken et al. 1974;McKay et al. 1978;Arndt et al. 1984;Delano 1986), combined with the high magma discharge rates (Wilson & Head 2017a), ensure that most of them do not cool significantly and coalesce in a lake feeding lava flows.Exposure of the pyroclasts to the lunar vacuum ensures extremely efficient release of volatiles, and the resulting lavas are expected to be almost completely vesicle free.We now explore the consequences of these conditions.

Conditions in the Lava Fountain
The height, H, of a lava fountain is given by the classical formula where U p is the clast eruption speed after the acceleration phase occupying a small region above the vent (Morgan et al. 2021) and g is the acceleration due to lunar gravity, 1.62 m s −2 .The maximum ballistic range of the pyroclasts, R, is given by implying that the height of the lava fountain is one quarter of its basal diameter.Note that Equations (2) and (3) neglect the curvature of the Moon's surface, a good approximation since R and H are always more than two orders of magnitude smaller than the lunar radius.The speed U p in the above equations is directly related to the mass fraction, n, of volatiles in the magma whose release and expansion into the vacuum above the surface are responsible for the dispersal of the pyroclasts.A good approximation (Wilson & Head 2017a) is where Q is the universal gas constant, 8.314 kJ kmol −1 K −1 , T is the magma temperature, ∼1600-1700 K for lunar basalts, and P f and P d are, respectively, the pressure at which the magma is disrupted into pyroclasts and the pressure at which pyroclasts decouple from the gas expansion.With lunar volatiles being dominated by CO and H 2 O in approximately equal proportions, the effective molecular weight, m, of the gas mixture is ∼0.5 (28 + 18) = 23 kg kmol −1 .More elaborate treatments taking account of thermodynamic reactions between the gas species released are possible; for example, the compositions suggested by Gaillard & Scaillet (2014) and Renggli et al. (2017) would imply values of m = 27 and 32, respectively, but there is considerable uncertainty as to which of these treatments might be more applicable.These issues are discussed by Wilson & Head (2017a), who showed that for a monodisperse distribution of pyroclast grain sizes the relationship between V and n was expected to be ∼U p = (5.278n) 1/2 .More recently, Morgan et al. (2021) modeled the relationship in more detail, finding that the clast size distribution should be bidisperse, giving rise to the relationship where U p is measured in meters per second and n in parts per million by mass.This implies that the maximum pyroclast range is given by where R is measured in meters, and again n is in parts per million by mass.Morgan et al. (2021) summarize published estimates of the total volatile content of lunar magmas; the largest estimate is 3400 ppm (Rutherford et al. 2017), which implies values of R up to 20 km.A point-source vent would be expected to produce a circularly symmetric fountain whose overall shape mimics the umbrella-like eruption "clouds," actually also lava fountains, seen on Io (Cataldo et al. 2002).Lunar examples would typically have been ∼100 times smaller than those on Io due to the lower volatile contents (except in rare special cases of secondary gas concentration in stalled dikes; Head et al. 2002).However, the shape of a fountain from a fissure vent depends on the ratio of the maximum range of pyroclasts, R, to the length along the strike of the fissure, L (Figure 4).If R is much greater than L (Figure 4(a)) the details of the shape of the vent will have a negligible effect on the distribution of clasts in the outer part of the fountain, and the number density of clasts will decrease as 1/R 2 in that region.In contrast, if R is much less than L (Figure 4(b)) then the number density of clasts will decrease more slowly as 1/R in all parts of the fountain except near the ends of the fissure.We saw above that R is likely to be as much as 20 km for the ∼3400 ppm total magma gas contents suggested by Rutherford et al. (2017) for a tenfold smaller gas content R would be as small as 2 km.We therefore expect the complete spectrum of possible eruption fountain geometries to have occurred on the Moon.
The number of clasts per unit volume in the fountain is directly proportional to the dense rock equivalent magma volume flux being erupted from the vent, F. The presence of large numbers of small droplets, with mean diameter d, in the lava fountain can lead to the clasts in the inner part of the fountain being prevented from radiating heat into space (their only potential method of cooling) due to the screening by clasts in the outer part of the fountain (Figure 2(c); see Head & Wilson 1989 for optical density effects in terrestrial Hawaiian eruptions).The thickness, X, of the shell that defines the outermost part of the fountain from which clasts can lose heat efficiently is shown by Wilson & Head (2017a) to be given by X d g R F 6.17 The numerical constants in these equations arise from the differing spatial distribution of clasts within the fountains.These equations allow us to calculate the fractions, f hot , of the pyroclasts that will fall to the surface from the fountain uncooled: for the point-source case and for the fissure-source case.
In general, there will be mixing, as clasts reach the ground, between those that have remained hot and those that have efficiently cooled.In modeling lunar pyroclast dispersal, Morgan et al. (2021) pointed out that the range of angles from the vertical at which clasts are ejected is much greater in eruptions in vacuum than in eruptions on Earth (e.g., Figure 1(d)).Morgan et al. (2021) suggested a maximum angle of 65°based on data related to supersonic gas jets; this is consistent with an analysis of pyroclastic deposits on Io by Glaze & Baloga (2000), who found that a value of 75°better fitted their data than the other two angles they modeled, 45°a nd 20°.These findings contrast strongly with, for example, our observations (Parfitt & Wilson 1994) of slightly pulsating vertical lava fountains during Phase 34 (August 1984) of the Pu'u'O'o eruption of Kilauea, Hawai'i, where pyroclasts were ejected at angles from the vertical ranging between 3°and 6°.These angles are much less than the ∼45°angle corresponding to the maximum range, and so pyroclasts in terrestrial fountains reach a given location on the ground having traveled by a path that began with a single launch angle.Furthermore, the clasts which land furthest from the vent are those that have spent the longest time in the outer part of the fountain from which heat is lost efficiently and are therefore the coolest, and commonly form cinder or spatter deposits (Head & Wilson 1989).In contrast, pyroclasts in fountains on the Moon can reach a given location on the ground via one of two paths.One path involves a launch angle less than the angle corresponding to the maximum range, which is a little less than 45°due to the dynamics of the clast-gas interaction (Morgan et al. 2021).The other path is more nearly horizontal and involves an angle greater than the critical value.The clasts in this second group spend much less time in the outer part of the fountain than those in the first group.The result is an intimate mixing, on the ground, of clasts that have very different thermal histories.We assume that the magma emerges from the vent at its liquidus temperature, T liq , and model the lava in the pond as having a uniform equilibrium temperature, T e, estimated by mixing the pyroclast fraction f hot erupted at the liquidus temperature T liq with the pyroclast fraction f cool = (1 − f hot ) that has cooled during flight to the ambient temperature, T a .The thermal equilibration takes account of both the sensible and latent heat exchanges as the cooled clasts are reheated to the solidus temperature, T sol , and then partially melted while the initially uncooled clasts partially crystallize as they cool, and leads to where Φ is the ratio L h /(T liq − T sol ).Here c c and c h are the specific heats of the lava in the subsolidus and melting ranges, respectively, and L h is the latent heat of fusion.In the same manner as Williams et al. (2000), we assume that the lava has the composition of the Apollo 12 low-titanium picritic basalt sample 12002, which has the property that only one species, olivine, grows as the lava cools from the liquidus temperature, 1713 K, to the solidus, 1423 K.The fraction of the lava that has crystallized, f c , therefore increases linearly from 0% to 100% as the temperature decreases over this range, so the fraction crystallized when the temperature reaches T e is given by Using Equation (11) for T e , Equation 12(a) becomes Substituting values for c c = 1000 J kg −1 K −1 , c h = 1570 J kg −1 K −1 , L h = 6.06 × 10 5 J kg −1 , and T a = 180 K for the Moon (Williams et al. 2017), assuming the solar heat flux ∼3 Ga ago was ∼20% less than now (Feulner 2012), T liq = 1713 K and T sol = 1423 K, we find that f c = 2.17 f cool .Other lava compositions would have differing thermal properties, but not by large amounts (Williams et al. 2000), and so it is a general result that cooling a given fraction of a mass of lava down to ambient temperature and mixing it back into the remaining hot lava causes close to double that fraction of crystallization after thermal equilibration.The fractions of pyroclasts landing uncooled, f hot , are shown as a function of the maximum pyroclast range in Figure 5 for a series of values of the erupted mass flux, F (curves labeled in cubic meters per second).Figure 5(a) is for a fissure vent with an outcrop length of 10 km, and Figure 5(b) is for a fissure of length of 3 km.For pyroclast ranges greater than ∼4 km, this shorter fissure acts as a point source rather than an elongate source.In each case, the upper axis gives the released magma gas fractions corresponding to the ranges on the lower axis.For the 10 km fissure, eruption rates are expected to be greater than 10 4 m 3 s −1 , implying that f hot is greater than 98%, so that f cool is less than 2% and, using Equation 12(b), less than f c = 4% crystallization will have occurred in the lava by the time it leaves the pond surrounding the vent.The situation is more complex if the fissure length has shrunk to 3 km as in Figure 5(b).In this case, if the released magma gas content is less than ∼800 ppm, so that pyroclasts land at up to ∼5 km from the vent, f hot is greater than 99%, so that f cool is less than 1% and no more than 2% crystallization occurs as lava passed through the fountain.However, at higher gas contents and clast ranges, progressively more of the clasts in the outer part of the fountain lose heat.For example, if F = 10 4 m 3 s −1 and the maximum clast range is 10 km, f hot decreases to 60%, f cool increases to 40%, and 80% crystallization can occur.An eruption with these parameters should form a substantial 10 km radius spatter rampart around a more nearly circular lava pond.Extreme examples of this may be the three cones at the south end of the vent structure of the sinuous rille Rima Schröter (Jawin et al. 2016).However, we stress that for the range of high-volume-flux conditions expected in Phase 2 eruptions, the temperature of lava leaving the vent pond will always be very close to its temperature as it approaches the surface from depth.

Implications for Lava Rheology
The quantity f c determines the initial rheology of the lava.Mafic magma erupted at its liquidus temperature and containing no gas bubbles can safely be assumed to be Newtonian, but even a small crystal content will induce non-Newtonian behavior (Dobran 2001).The simplest non-Newtonian fluid is a Bingham plastic, possessing a finite yield strength, Y, and a plastic viscosity, η p (Skelland 1967).Bingham rheology has been assumed by previous authors for flow modeling, but only for laminar flows (e.g., Dragoni et al. 1992;Dragoni & Tallarico 1994;Tallarico & Dragoni 2000).Here we explore its importance in turbulent flow.
For the functional variation of the liquid viscosity with temperature we use a power-law formulation as suggested by Hulme (1973): This function fits the viscosity data on Apollo sample 12002 given by Williams et al. (2000) with B = 1477.15K and q = 10.0207,making the liquidus viscosity 0.227 Pa s.We then evaluate the factor, f v , by which this value of η must be increased to give the plastic viscosity η p , due to the current value of f c .Williams et al. (2000) suggest combining the Roscoe-Einstein equation,  (2010).We adopt the following expression, which approximates the geometric mean of these data: with the values A = 0.04456 Pa, f o = 0.02209, f m = 0.6, and p = 1.77375.This expression implies yield strengths of a few hundred pascals at moderate solids contents, consistent with various measurements made on mafic lavas on Earth (Robert et al. 2014;Sehlke et al. 2014).Figure 6 shows how the yield strength Y and the viscosity correction factor f v vary with the lava crystal content f c ; also shown is the plastic viscosity η p of lava with the composition of the Apollo 12 basalt.23) shows that the yield strength is zero for crystal contents less than 2.1% and infinite for values greater than 45%, so the lines indicating the pond sizes are shown solid over this range of crystallinity.However, the viscosity is nonzero for all crystal contents and only becomes infinite when the crystal content is greater than 60%, so the pond size lines are shown broken for the highest and lowest crystal contents.
Various results follow from Figure 7. First, we expect the early Hawaiian phase of a lunar eruption to involve a fissure vent (Wilson & Head 2018;Figures 2(a), (c)), so Figure 7(a) is relevant.Morgan et al. (2021) show that, using current estimates of the volatile contents of lunar magmas, the maximum range of pyroclasts in lunar Hawaiian lava fountains is likely to be 20 km, so we take this as the maximum value of the pond half-width, R. Figure 7(a) shows that, for all volume fluxes greater than ∼3 × 10 4 m 3 s −1 , the lava in a 20 km halfwidth pond feeding initial lava flows will have a temperature very close to the liquidus, a crystal content less than the critical value of 2.1%, zero yield strength, and hence will exhibit Newtonian rheology.This is true for even smaller eruption rates if the pond width is smaller.However, we expect (Figures 2(a), (c)) that during the Hawaiian phase the volume flux will be much larger than 3 × 10 4 m 3 s −1 , decreasing through the range of values ∼10 6 to 10 5 m 3 s −1 .The Hawaiian phase is therefore guaranteed to produce near-liquidus-temperature Newtonian lava feeding the first flows leaving the vent pond.These flows consist of lava that has lost volatiles in the lava fountain and is expected to be essentially nonvesicular in nature.The density, ρ, of the Apollo 12 low-titanium basalt used in our modeling is 2900 kg m −3 , and Figure 7(a) shows that its viscosity, η, will be ∼0.2-0.3Pa s.We can therefore evaluate the Reynolds number of the lava, Re.For flows that are much wider than they are thick, true for all of the wellpreserved lunar flows in Mare Imbrium for which we can observe the three-dimensional morphology, the equivalent diameter of the flow is 4 times the flow thickness (Hulme 1973;Pinkerton & Wilson 1994) and so Re is given by where D f and U f are the flow thickness and mean speed, and we use η rather than η p because the lava will still be Newtonian at this point.Although at this stage we do not know the individual values of U f and D f , we do know their product, because continuity requires the equivalent of Equation (1), i.e., (U f D f ) = (F/L), and we have already shown that L will be ∼10 km.So, for values of F of 10 6 , 10 5 , 10 4 , and 10 3 m 3 s −1 , the values of Re will be close to 4 × 10 6 , 4 × 10 5 , 4 × 10 4 , and 4 × 10 3 , respectively.For Newtonian fluids, the motion is turbulent if Re is greater than ∼2200, and so, given the expectation that F will be much greater than 10 4 m 3 s −1 , it is guaranteed that the lava flows leaving the vent pond formed by a fissure vent will be turbulent.
The situation is less extreme for point-source vents in Figure 7(b).The spatial density of fountain pyroclasts traveling radially away from a point-source vent decreases as the cube of the distance traveled, and so much greater cooling can occur.For the highest magma volatile contents, ∼3000 ppm expelling pyroclasts to 20 km, an erupted volume flux greater than 10 6 m 3 s −1 would be needed to ensure negligible clast cooling, though much smaller fluxes would be needed for smaller volatile contents, e.g., 3 × 10 4 m 3 s −1 for 300 ppm.However, we would not expect an initially 10 km long fissure to have contracted to the extent of forming a near-circular point-source vent until the very late stages of an eruption.Although the volume flux would have decreased by then, it is likely that the activity would have become Strombolian (Section 3.1; Wilson & Head 2018) and Figure 7(b) would no longer be relevant.Hulme (1973) suggested that at least some lunar eruptions involved lava flowing in a turbulent manner to explain the formation of sinuous rilles.The above analysis implies that turbulence was near universal in lunar lava flows leaving the vent during the Hawaiian-style fissure-fed phase of lunar eruptions.

Morphology of Turbulent Flows
Although we can find no documented records of observations of an unambiguously fully turbulent lava flow on Earth, turbulent flows have been invoked to explain the thermal signatures of lava flows on Io (Williams et al. 2001;Davies et al. 2020).We consider that such flows must differ in fundamental ways from laminar flows whether they have Newtonian or non-Newtonian rheology (Skelland 1967).Figures 8(a The fact that shear stress increases with depth in an unconfined flow on an inclined plane makes it easy to understand why the presence of a yield strength causes the formation of a plug in the upper part of a laminar flow (Figure 8(b)).The expected internal structure of a turbulent Bingham flow is less clear, however.We infer that, during full turbulence, pockets of unsheared fluid rotate within intensely sheared regions; these unsheared pockets coalesce to form the plug as the motion transitions to laminar. Figure 8(e) gives a simplified impression of this process and shows the thermal boundary layers that begin to form at the top and base of the flow only after it becomes laminar.
Given the consequent continuous disturbance of the flow surface during turbulence, we support Hulme's (1973) argument that a strongly cooled crust cannot be formed on the surface of a fully turbulent lava flow.Support for this comes from thermal measurements on large-scale lava flows on Io (Davies et al. 2020).The extremely large peak thermal fluxes (>10 12 W μm −1 ) and rapid exponential decay times (days) measured for these eruptions can only be explained by radiation from initially turbulent lava flows.We therefore disagree with models of turbulent lava motion in komatiites (Huppert & Sparks 1985;Williams et al. 1998) and in lunar (Williams et al. 2000), Martian (Williams et al. 2005), and Ionian (Williams et al. 2001) lava flows which assume that a cooled crust can be stable on the surface of a turbulent flow.
Certainly, lava at the top surface of a flow begins to cool during the time it resides there, but it is efficiently mixed into the interior where it is reheated for a significantly longer time at the expense of the temperature of the bulk of the flow.A patch of lava on the exposed surface of an 8 m thick turbulent flow advancing at 4 m s −1 (values that we shall see in Section 3.3.9are plausible for lunar flows) will typically be exposed for a time, τ e , of ∼(8 m/4 m s −1 =) 2 s, in which time it will chill to a depth of ∼(2.3 [κ τ e ] 1/2 =) 3.25 mm, where κ is the thermal diffusivity, ∼10 −6 m 2 s −1 (Turcotte & Schubert 2002).The chilled platelet of lava will move down, then laterally, and finally upward in a convection cell in the flow interior for ∼3 times longer than it was exposed on the surface.During this time, it will have all its faces in contact with lava and thus will readily be reheated to the mean internal temperature.Thus, the entire bulk of the flow will cool very nearly uniformly provided full turbulence is maintained.
It is not immediately clear what limits the initial lateral spreading of turbulent flows of lava erupted at near-liquidus temperatures with minimal crystal contents and negligible yield strengths.Laminar flows on Earth (Figure 1(e)) deposit a layer of cooled crustal material onto surfaces over which they pass and incorporate very little surface material.Their levees consist of material that has been pushed aside to the edges of the flow by the advance of the flow front and, commonly, of material added later as the flow develops (Glaze et al. 2009).In general, such levees have a very different rheology from the uncooled lava in the core of the flow.In contrast, we argue above that turbulence prevents selective cooling of a turbulent flow surface, so that no cooled material is available to be deposited.Also, a turbulent flow has only a thin laminar boundary layer at its base and is more likely to incorporate loose fragments.This is particularly important on bodies without atmospheres, like the Moon, which develop fragmental regoliths on their surfaces due to unhindered meteoroid bombardment.The advancing boundaries of a turbulent lava flow may therefore pluck loose fine-grained regolith from the underlying surface and incorporate it into the body of the flow to add to the solids already present, modifying the rheology.Estimates of regolith formation rates during the main formation period of mare flows are ∼5 mm Ma −1 (Horz et al. 1991).Combined with typical intervals between eruptions in a given region, ∼20 Ma (Hiesinger et al. 2003), this implies that ∼0.1 m of regolith would commonly have been present.However, if the intervals between eruptions in the same location are short, very little regolith will have accumulated.Nevertheless, layers up to at least a few centimeters thick of unwelded fragmented glass shards may have formed at the tops of highly vesicular flows due to explosive auto-brecciation and so may be available for incorporation into turbulent flows (Head & Wilson 2020).
To explore this process, we assume that a flow spreads laterally due to vortices in the motion migrating sideways away from the downslope axis of the flow.If we assume that our typical 8 m thick flow spreads laterally and downslope at the same speed, 4 m s −1 , and that the largest turbulent vortices at the edge of the flow are cylindrical and equal in diameter to the thickness of the flow, then the perimeter and cross-sectional area of the vortex are (π × 8 =) 25.1 m and (0.25 π × 8 2 =) 50.3 m 2 , respectively.In one rotation, the vortex, and hence the edge of the flow, advances 25.1 m.If it accumulates 3% of a 0.1 m thick layer of clasts, i.e., a 3 mm thick layer, then these represent a cross-sectional area of (25.1 m × 0.003 m =) 0.075 4 m 2 , which is a fraction (0.0754/50.3=) 0.001 5, i.e., 0.15% of the vortex area and, since it is cylindrical, its volume.The added clasts, which we assume are derived from regolith formation on an earlier flow, are heated by the lava and melt, in turn cooling the bulk lava so that additional olivine crystals form.To get a conservative estimate of the importance of the process, we take an extreme case where the lava is initially crystal free at its liquidus T liq and entrains a volume fraction f r of regolith at the ambient temperature T a .Taking account of the exchanges of sensible and latent heats as in the case of Equation (11), the equilibrium temperature T e in this case is where, as before, Φ is equal to L h /(T liq − T sol ), and the fraction of crystals produced in the lava is given by Equation 12(a), i.e., f c = (T liq − T e ) / (T liq − T sol ).The final solid content of the resulting mixture is the sum of the lava that has crystallized and the unmelted fraction of the added regolith, and the final liquid content is the sum of the melted regolith and uncrystallized lava.The solid fraction of the total is denoted f s .When lava incorporates regolith in this way, it is the total solids content, f s , that controls the rheology, so f s replaces f c in Equations 14(a) and 15(a).Table 1 shows how T e decreases and f s increases with f r .f s reaches 45%, which Equation 15(b) shows is sufficient to produce an infinite yield strength, when f r has increased to 28.6%.In practice, f s = 43% would produce a yield strength of 75 Pa, sufficient to halt an 8 m thick flow on a sin α = 0.002 (∼0.1 degree) slope, requiring f r = ∼27%, a condition that would require (27/0.15=) 180 vortex rotations and be reached after the edge of the flow had spread (180 × 25.1 =) 4518 m, i.e., ∼4.5 km.Thus, taking account of both levees, our nominal 10 km wide flow would have spread to a total width of [10 + (2 × 4.5) =] ∼19 km after its front had traveled 4.5 km downslope and would then have ceased to spread laterally.If now we assume that 3 cm of fragmental material is assimilated at each vortex rotation instead of 3 mm, this would represent 1.5% of the vortex volume, and only (27/1.5 =) 18 vortex rotations would be needed to reach f r = ∼27%.Each flow edge would have spread by (18 × 25.1 =) 452 m, producing a 10.9 km wide flow that ceased to spread laterally after its front had traveled 452 m downslope.We could plausibly assume other values for the initial ratio of lateral spreading rate to downslope advance speed; indeed, the ratio likely depends on both the ground slope and the flow thickness.Nevertheless, even if the magma forming a turbulent lunar lava flow is erupted crystal free, ingestion of a small amount of regolith material could quickly induce a large enough yield strength to form a levee able to limit the lateral spread of the flow.
A second factor that may contribute to the generation of a yield strength in lava at the edge of a flow near the vent is controlled by the cooling of pyroclasts in the fire fountain.Section 3.3.1 pointed out that, depending on their launch angle, different pyroclasts will spend different fractions of their trajectories deep within the opaque, thermally isolated part of the fountain or in the zone near the edge where they will cool.The pyroclasts launched at angles close to the angle giving the maximum range will spend a disproportionate length of time in the cooling zone and will be the coolest of all those reaching the ground.Whatever the detailed geometry of lava leaving the vent pond to form a flow, it is this material that will be incorporated into the edges of the flow, where it will contribute to increasing the initial yield strength of the marginal lava forming levees.At present it is difficult to quantify this pondgenerated contribution.

Effect of Rheology on Flow Dynamics
Lava cools as it travels away from the vent.We model the detailed thermal history in Section 3.3.5,but first use the discussion of rheology in Section 3.3.2 to explore whether a flow is likely to continue to be turbulent as its crystal content increases and it becomes non-Newtonian.If the erupted volume flux remains large and the viscosity remains small, the flow will be turbulent.However, the volume flux is expected to decrease as the eruption progresses (Wilson & Head 2017a, 2018; Figure 2(a)) and the viscosity will certainly increase as cooling increases the crystal content.It is thus not possible to establish a priori if the flow will be turbulent or laminar.Wilson & Head (1981) pointed out that, for laminar flows, if the flow speed is found from the formulae for both turbulent and laminar flow, then whichever formula gives the smaller speed will always be the correct one to use.
We now consider this issue for Bingham plastic fluids.If the flow is fully turbulent, the speed, U t , is given by where D t is the flow thickness, α is the slope of the underlying ground, and the Fanning friction factor, f f , can be expressed in terms of the flow Reynolds number Re by (Skelland 1967) f R e 0.0014 0.125 , 19 where Re is now defined in terms of the plastic viscosity η p by The initial width of the flow is likely to be similar to the length of the fissure, so continuity requires that which gives us a second relationship between U t and D t .Equations ( 18)-( 21) can be combined to eliminate U t and give D t analytically in terms of F, L, ρ, α, and η p :

= ( )
This relationship is valid for He > ∼5000; at smaller values the rheology is extremely close to Newtonian.If the lava is erupted at a temperature significantly below its liquidus, or if the eruption rate is small enough that a significant fraction of the clasts in the lava fountain are able to cool in flight, the lava leaving the vent pond to form a flow will have a finite crystal content, a nonzero yield strength, and a viscosity greater than the liquidus value.
Even if the lava leaves the vent at its liquidus it will cool and crystallize as it flows away.There will therefore be finite values for the Hedström number He and the critical Reynolds number Re crit , while at the same time the higher viscosity causes the flow speed and the actual Reynolds number Re to be smaller than if the lava were hotter.Nevertheless, provided Re remains greater than Re crit , the motion is turbulent.
If Re becomes less than Re crit , however, the motion becomes laminar.Cooling at the flow margins will likely maintain the laminar state unless some other factor, for example a large increase in the slope of the ground beneath the flow (Glaze et al. 2014), intervenes.The main characteristic of a Bingham plastic fluid flowing in a laminar fashion is the presence of an unsheared layer of fluid, commonly described as a plug, extending from the upper surface down to the level where the shear stress is equal to the yield strength (Figure 8(e)).The thickness of this layer, h, is given by and we define the ratio, j, of this plug thickness to the total thickness of the laminar flow, D l : The mean lava speed in laminar flow is then given by (Skelland 1967) as If the yield strength Y is zero there is no unsheared plug, so h is zero, j is zero, and Equation and U l follows from Equation (21).The general case of nonzero yield strength and plastic viscosity is more complicated, however.By combining Equations ( 27), (26), and (25), we can obtain a cubic equation in D l : where Ω = Y/(ρ g sin α).Various analytical solutions of cubic equations exist but all are computationally tedious.A simple solution is to rewrite Equation 29(a) as If we take an initial estimate of D l and insert this into the righthand side of this equation, the resulting new value of D l generated is a better estimate of the true value; the process is then repeated.The thickness of a flow must always be greater than the thickness of its unsheared plug, and so whenever the yield strength is nonzero, an initial estimate equal to the plug thickness given by Equation (25) yields rapid convergence, typically better than 1 part in 10 6 in six iterations.When Y is exactly zero, a value of 10 m is used as the initial trial for D l , on the grounds that the observational evidence presented in Section 2 suggests that mare lava flow thicknesses are expected to lie in the range 10-20 m.We now have two criteria for judging when the transition from turbulent to laminar flow occurs in Bingham plastic lava.We either require Re from Equation (20) to become equal to Re crit from Equation (24) or, as in the case of Newtonian fluids, require U t from Equation (18) to become equal to U l from Equation (27).The two criteria should result in the same transition conditions but, because they rely on fitting functions to experimental data only available as graphs, they cannot be expected to agree perfectly.In practice, we find that the values of the distance from the vent at which the turbulent-to-laminar transition occurs calculated from the two approaches differ by no more than 2%.An additional error is caused by approximating experimental measurements of f f as a function of Re and He by two separate functions at low and high values of Re.The data plotted in Figures 5.5 and 6.10 in Skelland (1967) imply that the theoretical curves in his Figure 3.2 must always slightly underestimate the friction factor over a range of Reynolds numbers spanning the transition from laminar to turbulent motion.However, fitting a smooth curve by eye through a sample of the data points and comparing the implied flow speeds with those from the fitted functions shows that the error in the distance from the vent of the transition point is less than 1%.
3.3.6.Thermal Evolution of the Flowing Lava Hulme (1973) showed that heat loss from a turbulent lava flow would be dominated by radiation from the upper surface of the flow, continuously replaced by material stirred from below by the turbulence.The radiated energy loss per unit mass, ΔE R , as the flow advances a distance Δx is given by where s is the Stefan-Boltzmann constant, 5.67 × 10 −8 W m −2 K −4 , and ε, U f , and D f are the emissivity, mean speed, and thickness of the lava, respectively.We assume ε = 0.96 (Davies et al. 2011) for the high temperatures involved.Hulme's (1973) original work was aimed at explaining how turbulent flows were able to heat the ground beneath them to its melting point, thus explaining the eroded channels of lunar sinuous rilles.Essentially the same models were expanded by Huppert & Sparks (1985) and Williams et al. (1998) to treat channels eroded by komatiite lavas on Earth and by Williams et al. (2000) to revisit lunar rilles.All these treatments calculated heat transfer to the ground on the assumption that melting was already taking place.However, Hulme (1982) pointed out that a finite time was needed for the ground to heat to the point of melting.Assuming the lava is flowing over earlier flows of similar composition, the melting temperature is the solidus temperature T sol .Hulme (1982) shows that the time needed to initiate melting at any given distance from the vent, τ m , can be related to the time that lava takes to flow from the vent to this location, τ f , by where Λ is the thickness of the heated layer in the ground beneath the flow divided by the thickness of the thermal boundary layer at the base of the flow, and is given by Hulme (1982) as where T and T sol are the bulk lava temperature and the solidus as before, and T a is the average ambient temperature, again ∼180 K. Defining T r as the ratio (T sol − T a )/(T − T a ), this is more conveniently reorganized as To give quantitative examples of the importance of this issue we again assume the thermal parameters for the Apollo 12 sample 12002.For lava having just left the vent at a temperature close to the liquidus, say 1700 K, Λ is 6.733, τ m /τ f is 28.9, and so at a location 100 m from the vent, reached after 20 s by lava flowing at 5 m s −1 , melting starts after only 9.6 minutes.Anticipating the results from the flow models described below in Section 3.3.9,if the eruption rate is near the upper end of the range inferred by Wilson & Head (2017a) for major mare lava flow eruptions, 10 6 m 3 s −1 , and decreases exponentially with a decay constant of 24 hr, the transition from turbulent to laminar flow occurs at 132.6 km from the vent 8.67 hr after the eruption begins.With a flux decay rate of 24 hr, the volume of lava erupted after 3 days is 82.1 km 3 , 95% of the total 86.4 km 3 of lava erupted.Thus, the bulk of the lava was erupted before substrate melting and thermal erosion could have started in all but the most proximal parts of the flow.Other than very near the vent, a negligible amount of the flow's heat was used to melt its substrate.In that case the heat transfer to the ground per unit mass of lava is adequately modeled by (Hulme 1973) where U f and D f are the mean speed and thickness of the lava flow irrespective of whether the motion is turbulent or laminar.
Here h C is a heat transfer coefficient, given by where k is the thermal conductivity of the lava, and Pr is the Prandtl number for the motion, given by Pr .36 The energy per unit mass produced by viscous dissipation as work is done against friction with the ground is where α is the slope of the ground over which the flow is moving.
Finally, the energy per unit mass released as the temperature decreases by ΔT during crystallization is where c h and L h are the specific heat and latent heat of the lava, respectively.Equating the heat generated to the heat lost, where again Φ = [L h /(T liq − T sol )], which allows the temperature evolution to be followed as the flow advances.As the temperature decreases, Equations (13)-( 15a) are used to evaluate the crystal content and hence the yield strength and plastic viscosity, and then Equations ( 18)-( 29a) are used to find the flow thickness and mean speed, taking account of the turbulent-to-laminar transition when this occurs.

Limitations on Flow Front Advance
While a lava flow remains turbulent, we approximate the temperature as being uniform throughout the flow.However, when turbulence ceases, thermal boundary layers grow at the top and base of the flow.We approximate the thickness, z, of these layers by the solution to the Stefan problem (Turcotte & Schubert 2002) so that where τ c is the time since this cooling began at the turbulent-tolaminar transition.The boundary layer at the top of the flow initially grows downward into the molten lava plug, of thickness h given by Equation (25), which forms the upper part of the flow.The plug is already incapable of being sheared due to the existence of the yield strength.As a result, its solidification as the wave of cooling passes through it has no influence on the flow's motion until the thickness of the boundary layer exceeds the plug thickness.In practice, we find that this never happens for any of the modeled flows, and that z is always very much less than h.At the base of the flow, which is still being sheared, the boundary layer growth plays the same role as in a fully laminar flow.Thus, we assume that, as the basal cooling wave propagates into the flow, the additional crystals that form in the already partly crystallized flow adhere to the substrate and are lost from the moving part of the flow., or some or all parts of the flow deposit inflate as lava is injected into the flow's interior, or both processes occur.This combination of rheological and thermal limitations on the maximum travel distance of a lava flow involving a turbulent phase has various consequences.In a fully laminar flow it is the meeting, near the middle of the flow, of the propagating upper and lower thermal fronts that limits the advance of the flow (Pinkerton & Wilson 1994).However, we find that in our initially turbulent modeled lunar flows, the thickness of the rheologically controlled unsheared plug that forms after turbulence ceases is always much greater than the thickness of the upper thermal boundary layer due to conduction (Figure 8(e)).As a result, the thickness of the layer of lava that can be sheared is much smaller than if it were controlled by the convergence of two purely thermal boundary layers.This means that an initially turbulent flow will not travel as far as a laminar flow that is compositionally similar and erupted at the same volume flux per unit flow width.Nevertheless, in general, flows that exhibit turbulence do so because their lava has a low viscosity and is erupted at a high volume flux, and this combination of factors allows such flows to reach lengths similar to or even greater than those reached by laminar flows despite the fact that the turbulent flows have high cooling rates while turbulence persists.
A second predictable consequence of the high initial cooling rates of turbulent flows is that such flows should contain two classes of crystals.We shall see in Section 3.3.9that turbulence commonly stops when the crystal content of the lava has grown to ∼40%.The crystals that grow during the turbulent phase do so with cooling rates of ∼10-20 K hr −1 .The remaining 60% of the crystallization occurs at much slower rates, the value varying with position in the flow and the total flow thickness.
Finally, the fact that high values of crystallinity, viscosity, and yield strength are all generated quickly during the turbulent phase of a flow minimizes the opportunity for fractional crystallization to occur.Formation of small gas bubbles by second boiling later in the cooling and crystallization process further inhibits crystal settling.

Implementation of Model
The above system of equations has been organized as a spreadsheet.The temperature of the lava leaving the vent is taken to be 1700 K, a little less than the liquidus temperature, 1713 K.This is to allow for the possible assimilation of small amounts of dike wall rock during magma ascent, and inclusion of preexisting regolith formed by impacts or by autobrecciation (spontaneous disintegration under thermal stresses; Head & Wilson 2020) of previous vesicular lava flows.The thickness of the regolith layer incorporated into the flow front is assumed to be 5 cm.The initial conditions for a given model run are the initial total volume flux of lava from the vent, F i , the time constant for the exponential decrease of the volume flux with time, τ d , the length of the vent along the strike, L, which also defines the initial width of the lava flow leaving the vent, and the slope, α, of the ground across which the lava flows away from the vent.
In modeling the supply of magma from mantle sources to the lunar surface, Wilson & Head (2017a) predicted that largevolume eruptions in the maria should involve an initial magma volume flux of a few × 10 6 m 3 s −1 , the flux decaying on a timescale of a few days in a manner that depended in detail on the rheological response of the lithosphere to the arrival and subsequent horizontal closure of the dike carrying the magma.Early models of sinuous rille-forming eruptions (Head & Wilson 1981) implied initial fluxes of a few × 10 4 m 3 s −1 , with much longer decay constants (weeks to months) being needed to allow time for extensive thermal erosion.Since we are not concerned here with eruptions causing appreciable erosion, we use initial lava volume fluxes ranging from 3 × 10 6 to 10 4 m 3 s −1 , and for mathematical simplicity we assume an exponential decay of the flux, with e-folding times ranging from 0.5 to 4 days.With this assumption, the volume, V, of lava erupted after a time t is and as t goes to infinity the total erupted volume, V tot , is just The Wilson & Head (2017a) model of lunar eruptions links erupted volume flux to dike volume; given the common ∼10 −3 aspect ratio of dikes, this suggests a link between the initial flux F i and the initial fissure vent length L. We therefore explore various permutations of the above fluxes with fissure lengths ranging from 0.5 km for the lowest flux eruptions to 40 km for the highest.
The baseline ground slope used is 0.002 radians, close to the typical values in the southwest part of Mare Imbrium hosting its well-preserved late-stage flows (Chen et al. 2021).However, this may not have been the slope at the time of the eruptions.Loading may have caused the center of the mare basin to sink (Solomon & Head 1979, 1980), or isostatic mantle recovery from the original impact may have caused it to rise (e.g., Freed et al. 2014).We therefore run models with slopes ranging from 0.001 to 0.004 radians (∼0.06°to ∼0.2°).
For various permutations of the above initial conditions, we follow the motion of a batch of lava from the vent to the flow front as the flow advances and the supply rate at the vent declines.The parameters evaluated as a function of time and distance are speed, thickness, temperature, crystal content, total solids content, yield strength, viscosity, Reynolds number, Hedström number, and critical Reynolds number.After the transition from turbulent to laminar flow, we evaluate the thickness of the unsheared plug forming the upper part of the flow, the thickness of the cooling thermal boundary layer at the base of the flow, and the fraction of the thickness of the flow that still consists of lava that can be sheared and so allow flow advance.After the flow front comes to rest, the difference between the lava volume erupted so far and the amount of lava still to be erupted is used to find the percentage by which the stationary flow must inflate, using the assumption that postinflation second boiling (e.g., Edmonds & Woods 2018) as the flow cools (Wilson et al. 2019) leads to the formation of a vesicle content of ∼75%.The results of 17 flow simulations are summarized in Table 2.

Discussion of Phase 2 Results
The four stages of mare basalt lava flow eruptions are shown in Figure 2(a), and the first of the effusive stages, Phase 2, is shown in Figure 2(c).The greatest flow lengths predicted, in the range 500-750 km, are similar to the largest estimates for the well-preserved flows in Mare Imbrium (Schaber 1973), and the thicknesses, up to at least 13 m, encompass the ∼3 to 14 m range estimated from exposures of flows in craters walls (Robinson et al. 2012;Rumpf et al. 2020).Figure 9 shows the critical variables for flow model #1 in Table 2, an eruption from a 20 km long fissure vent onto a slope of 0.002 radians with an initial eruption rate of 10 6 m 3 s −1 and an exponential decay constant of 1 day.Lava leaves the vent as a 10 m thick, turbulent flow (Figure 9(b)) with a mean speed of 5 m s −1 (Figure 9(a)).It becomes laminar at 133 km from the vent, where the thickness is still close to 10 m and the speed has decreased to ∼3.4 m s −1 , and reaches a maximum length of 421 km. Figure 9(c) shows the temperature of the well-mixed turbulent flow decreasing from a near-liquidus 1700 to 1592 K as it becomes laminar.The decreasing temperature increases the solids content (Figure 9(d)), the yield strength (Figure 9(e)), and the viscosity (Figure 9(f)).The core of the flow remains at 1592 K, and its solids content and rheology remain fixed, as thermal boundary layers at the top and base of the flow become thicker (Figure 9(g)).Figure 9(h) shows the vertical extent of the flow that cannot be sheared because of the presence of the yield strength.In the laminar part of the flow, the region more than 133 km from the vent, this unsheared region is simply the plug forming the upper part of the flow, as shown in Figure 8(b).However, in the proximal, turbulent part of the flow this vertical extent does not represent a single discrete layer but instead must represent the sum of the vertical extents of a series of clots of unsheared fluid distributed randomly throughout the flow and tumbling within the general turbulence.Figure 9(i) shows the thickness of the layer of lava that can be sheared, found by subtracting the unsheared thickness components in Figures 9(g) and (h) from the total flow thickness in Figure 9(b).Flow advance stops when the thickness of this deformable layer becomes zero.Finally, Figure 9(j) shows how the rheological parameters, expressed through the Reynolds and Hedström numbers, control the flow dynamics.The Hedström number is initially ∼3 × 10 6 and first increases and then decreases as the the yield strength and viscosity of the lava evolve.The critical Reynolds number quite closely tracks the Hedström number on the logarithmic scale of Figure 9(j) due to the approximately square root dependence seen in Equation (24) and is initially ∼4 × 10 4 .The Reynolds number is initially ∼1.8 × 10 6 , much greater than the critical Reynolds number, thus guaranteeing turbulence, and decreases steadily until it becomes equal to the critical Reynolds number, at which point turbulence ceases.Thereafter the Reynolds number is always less than the critical Reynolds number, though initially only slightly.The model results given here are for lava moving down a uniform slope.Real planetary topography is irregular, and suitable variations of slope could cause both the Reynolds number and its critical counterpart to increase or decrease slightly while maintaining nearly the same values for a significant distance, opening the possibility of a temporary change back to turbulence.Changes of this kind caused by sudden variations of substrate slope or flow direction were suggested by Glaze et al. (2014).Table 2 demonstrates that the maximum lengths and typical thicknesses of flows vary widely depending on the initial volume eruption rate F, the flux decay constant τ d , and the fissure length L. We now identify the most important interdependencies.
Figure 10 shows how the maximum length of a flow, Z max , and the distance from the vent at which the turbulent-to-laminar transition takes place, Z trans , vary with the ratio F/L (Figure 10 2, that eruptions with numerous combinations of these three variables could have produced the up to ∼400 km lengths (Schaber 1973;Chen et al. 2018) of the well-preserved flows in Mare Imbrium.In contrast, the distance from the vent at which the turbulent-to-laminar transition takes place depends only weakly on τ d and α, and depends almost entirely on the value of F/L.This result will have particular significance for eruptions that continue long enough to cause thermal erosion of the substrate and sinuous rille formation.Since the width of the rille is a good proxy for L, combining this   with the length of the eroded channel should provide a good indicator of the eruption rate F.
Figure 11 shows how the total solids content, f trans , in the lava when the transition from turbulent to laminar flow takes place depends on F/L (Figure 11(a)), τ d (Figure 11(b)), and α (Figure 11(c)).The strongest dependency is on F/L and there is a significant contribution from α. f trans controls the plastic viscosity and yield strength of the lava in the laminar part of the flow, and these in turn control the thickness of the lava flow when its front comes to rest.F/L, τ d , and α also control the time at which the flow stops advancing.The sooner this occurs, the greater is the amount of lava still to be erupted, and it is this lava that leads to flow inflation.The net result of these factors can be illustrated in Figure 12, which shows the factor, I, by which the depth of the flow increases due to inflation after the flow front halts.The long curve in the figure shows how I varies with F/L for τ d = 1 day and α = 0.002.The two short segments give an indication of the effect of increasing τ d (and hence of increasing the total volume of lava erupted) to first 2 and then 3 days.The influence of changing α between 0.001 and 0.004 is very small and is indicated by the just-visible cluster of points above and below the curve at F/L = 30 m 2 s −1 .Note that the inflation shown in Figure 12 is that which would be caused by the injection of degassed, vesicle-free lava into a recently halted distal flow.Section 4 below describes additional inflation due to the consequences of the cooling of the intruded lava in cases where it has not already lost all of its volatiles.
Figure 13 shows the cooling rate, C t , of the lava in the turbulent part of the flow as a function of the ratio F/L.The cooling rate decreases slowly between the vent and the point of transition to laminar flow, typically by 16% ± 3%, and Figure 13 shows the average rate during turbulence for the models in Table 2.The scatter of data points for F/L = 30 represents the range for the data in Table 2, i.e., α between 0.0005 and 0.004, and τ d between 0.5 and 3 days.Cooling rates range from ∼10 K hr −1 for large values of F/L to several tens of kelvins per hour for small F/L.Small values of F/L arise  when either F is small or L is large.However, the eruption model of Wilson & Head (2017a) implies a positive correlation between F and L because they are linked via the feeding dike geometry.Combining Equations ( 5) and (8) in Wilson & Head (2017a), which give the functional dependencies of dike width and magma rise speed on dike spatial scale, with the continuity condition, our Equation (1), shows that L is proportional to F 2/5 .The highest eruption rate suggested by the Wilson & Head (2017a) model is ∼10 6 m 3 s −1 and the smallest is ∼10 4 m 3 s −1 .Raising the ratio of these eruption rates, 100, to the power 2/5 implies that the ratio of the longest to the shortest active fissures should be 6.3.The longest lunar fissure vent recognized by Head & Wilson (2017) extended for ∼13 km, implying that the shortest should be ∼2 km long.Recognizing such features unambiguously is difficult, but we note that in cases where the source ponds of sinuous rilles are elongate in shape (Oberbeck et al. 1971;Head & Wilson 1980, 1981), the difference between the long and short axes is of this order.These pairs of values of F and L imply that as L ranges from 2 to 13 km, F/L ranges from 5 to 77 m 2 s −1 .Figure 13 then shows that, during turbulence, the expected cooling rates range from 10 K hr −1 for high-flux eruptions from long fissures to as much as 40 K hr −1 for low-flux eruptions from short fissures.
After the motion becomes laminar, the lava in the core of the flow is insulated by the slowly growing thermal boundary layers at the top, base, and sides of the flow, and no further cooling occurs at any given depth in the flow until the thermal boundary layer reaches that depth.This stark contrast between the rapid cooling of the entire volume of the flow while it is turbulent followed by the protection of the core of the flow by its growing thermal boundary layers after it becomes laminar leads to an important result.The crystals in a flow that has experienced a turbulent phase should all show evidence for having initially cooled at the same high rate of 10-40 K hr −1 .The subsequent cooling rate recorded in the internal crystal growth pattern should be very much smaller and should vary widely between samples from different depths within the flow.Using a model of the cooling of basalt flows by Davies (1996), Table 3 gives typical initial cooling rates, C e , at various depths below the surface of a flow after it has been emplaced for a range of depths likely to be relevant for the Moon.Values are extremely large at the surface but at depths greater than 0.1 m they are always less than ∼10 K hr −1 , and at depths greater than 2 m are less than 0.04 K hr −1 , three orders of magnitude less than during the turbulent phase.

Phases 3 and 4
As the change from Hawaiian to Strombolian activity that constitutes Phase 3 occurs (Figures 2(a), (c), (d)), several changes take place.The volume eruption rate of magma decreases, and the width of the dike does not change greatly.This results in a decrease in the rise speed of the magma in the dike, allowing time for bubbles of CO nucleated in the deepest parts of the dike to expand, accelerate relative to the magmatic liquid, and overtake smaller CO bubbles nucleating at shallower depths.Coalescence accelerates the growing bubbles further and large gas bubbles arrive intermittently at the level of the lava lake around the vent.The magma through which they are rising is itself rising and, at depths less than ∼500 m (Rutherford et al. 2017), nucleating bubbles of H 2 O and sulfur compounds.It is inevitable that the large CO bubbles scavenge some of these other bubbles, so that the gas mass fraction in the bulk of the magma driving the fire fountain decreases.Also, the decompression of the CO bubbles creates large pockets of expanding gas that pass through the rest of the fountain.These pockets of gas do not couple strongly to the pyroclasts being created as the bubbles of H 2 O and sulfur compounds become   2, i.e., α between 0.0005 and 0.004, and τ d between 0.5 and 3 days.close-packed and the magma undergoes fragmentation, and so do not contribute much energy to the acceleration of these clasts.The net effect is to make the fountain become pulsatory and, because the total energy provided to the pyroclasts by gas expansion has decreased, to progressively reduce the height and width of the fountain.The change from Phase 3 to Phase 4 occurs as the lower tip of the dike feeding the eruption decelerates to zero as the dike reaches buoyancy equilibrium with the density structure of the crust and upper mantle (Wilson & Head 2017a).The magma rise speed by now has become very small, because the force driving magma upward is no longer the bulk motion of the rising dike itself, but rather the horizontal compression of the dike as the lithosphere dilation caused by the arrival of the dike relaxes.The very low magma rise speed greatly facilitates the overtaking of small gas bubbles by large ones and maximizes bubble coalescence.Eventually the largest bubbles will almost completely fill the width of the dike (Figure 2(d)), leaving only a thin film of magma against the walls, and can be regarded as gas slugs.The eruption style becomes that of Strombolian explosions punctuating the eruption of now very volatile-poor magma.The extent of this effect can be quantified using the results of Llewellin et al. (2012, their Figure 4), who treated fluids flowing in subvertical pipes where the scale is large enough that surface tension can be neglected, certainly the case for volcanic conduit and vent systems.They showed that the dimensionless wall film thickness, λ, defined by where l is the thickness of the liquid film against the pipe wall and r is the pipe radius, is a simple function of the dimensionless inverse viscosity of the liquid, N f , defined by where D p is the pipe diameter, ρ is the liquid density, and η is the liquid viscosity.For N f < 10, λ is independent of N f and has the value λ ≈ 0.33; in the interval 10 < N f < 10 4 , λ decreases sigmoidally with increasing N f ; and for N f > 10 4 , λ is again independent of N f and is ∼0.08.This treatment was developed for circular pipes rather than elongate volcanic fissures, but as long as λ is small, which we shall see is true here, the wall curvature can be neglected and D p can be equated to the dike width, W, which is equal to 30 to 100 m initially, narrowing to say 10-30 m (Wilson & Head 2017a) late in the eruption.The density of lunar magma can be taken as ∼2900 kg m −3 and the viscosity at the liquidus as ∼0.5 Pa s.This means that N f is initially ∼6 × 10 5 to 4 × 10 6 , reducing to ∼1 × 10 5 to 6 × 10 5 ; all of these values are much greater than 10 4 , so that λ is always ∼0.08, thus justifying equating D p to W. Then for D p = 10 to 30 m, i.e., r = 5 to 15 m, the film thickness l is 0.4 to 1.2 m.A film is present on each side of the dike, so the fraction of the cross-sectional area of the magma that interacts with the slug is (1-2 λ) = 0.84; we assume that all of the gas bubbles that are not preserved in the magma films on the walls are incorporated into the slugs, leaving the magma that reaches the surface with only 0.16, i.e., 16%, of its original volatiles.Rutherford et al. (2017) found that picritic lunar magmas could have released between ∼1100 and ∼2700 ppm total volatiles during their ascent; the upper value would decrease to 1500 ppm if most of the CO was lost during the Hawaiian phase of the eruption.Depending on the nature of the petrogenesis of the nonpicritic mare basalts, their water contents could be as much as a factor of 2 greater than in the picrites (M.J. Rutherford 2017, personal communication).Taken together, these results suggest that the maximum potential equivalent H 2 O content in lunar magmas could be ∼3000 ppm, and after scavenging by Strombolian slugs the residual amount would be ∼16% of this, i.e., ∼500 ppm, with a more likely value perhaps half of this.
No data exist on the solubility as a function of pressure of SO 2 and H 2 S in lunar magmas, but measurements exist for the solubility of H 2 O in some lunar-like magma compositions (Newcombe et al. 2017), and we use water as a proxy for the total volatiles in order to calculate the size and number of H 2 O bubbles formed in the last stages of magma ascent from the conduit into the lava lake filling the vent and feeding flows.Let the solubility as a function of pressure, P, be S(P).The data of Newcombe et al. (2017) where S is in parts per million and P is in pascals.If the total water content is n W , the mass fraction exsolved, n, is as long as this difference is positive and zero otherwise.Approximating H 2 O as an ideal gas, the specific volume fractions of released gas, v g , and liquid magma, v m , are given by where the magma temperature, T, is taken as 1500 K, and m is the molecular mass of H 2 O, 18.015 kg kmol −1 .The gas volume fraction, f g , in the bulk liquid-gas mixture, i.e., the vesicularity, is The rise speed of the magma in the dike by this stage of the eruption (Phase 3) is very small (Figure 2(a)), requiring only a very small pressure gradient to overcome wall friction, and so the pressure within the magma at any depth is very close to the static weight of the overlying magma column.This implies that the pressure, P, varies with depth y as P gdy, 4 9 where β is the bulk density of the magma, given by Equations ( 44) through (50) cannot be solved analytically but are readily implemented in a spreadsheet.They predict that the magma forms a foam with a vesicularity that increases toward the surface, where it is theoretically 100% because of the zero pressure.This impossible condition is avoided in practice because of surface tension.The pressure inside bubbles in any liquid is greater than the pressure in the surrounding liquid by an amount determined by the surface tension of the liquid-gas interface and the bubble size.Gas bubbles in magmas nucleate with diameters, f, of ∼10 μm at a depth determined by the mass fraction of the volatile, n, still remaining in the magma and the solubility law.Table 4(a) shows the diameter, f f , to which these bubbles of H 2 O gas will have expanded as they approach the surface of the rising magma in the lava lake.For n in the range 100-750 ppm the bubble diameters will have grown to diameters between ∼300 and 750 μm by the time they burst.The surface tension of the basalt-water interface, σ, is ∼0.37 J m −2 (Mangan & Cashman 1996), and so the excess internal pressure in bursting bubbles will be P s = 2σ/f, which Table 4(a) shows will range from 5.3 to 2.0 kPa.As long as a continuous liquid film surrounds each bubble, the bubble will be stable.
However, if the liquid film drains or becomes extremely thin due to close packing of bubbles, film collapse and interconnection of bubbles become likely.When this happens very close to the top of the magmatic foam the difference between the internal gas pressure and the vacuum will cause the top layer of bubbles to explode.The speeds, U p , to which the resulting magma droplets can be accelerated are again given by Equation (4), now with P d = P s and P f = 90 Pa as before.Table 4(a) shows that maximum speeds range from ∼16 to ∼32 m s −1 and that liquid droplets could be thrown to heights, H s , of 84-307 m.However, as lava moves radially outward from the center of the lake, interference between the products of the early-exploding layers descending and the products of the later-exploding layers rising will produce a much lower cloud of droplets.This problem was treated for explosive eruptions producing pyroclastic flows on Earth and Mars by Wilson & Heslop (1990).The mechanism is the same here, though with much smaller volatile amounts.We have implemented Equations (1)-( 4) of Wilson & Heslop (1990) as a spreadsheet, and Table 4(b) gives the resulting values of U p , H s , and the total time of flight of pyroclasts, τ str .We show below that lava will typically move away from the vent at speeds of ∼1 m s −1 during the Strombolian phase of the eruption, so the flight times of ∼10-30 s imply that stable porosity profiles will be attained after the lava has traveled ∼15-30 m.
The process described by Table 4 is best thought of as the rapid growth and subsequent slower decay of an expanded fluidized bed (Figure 14) as a given batch of lava moves away from the vent.The profile of the lava flow after the process is complete is expected to consist of an unwelded layer of shattered and chilled bubble wall fragments grading downward into a welded layer of such fragments, grading quite abruptly into very vesicular foam (Wilson & Head 2017b;Wilson et al. 2019, their Figure 5).The porosity of the foam decreases with increasing depth until it becomes zero at the depth where the residual volatile content of the magma is just saturated at the pressure corresponding to the weight of the overlying material.We assume that the boundary between stable foam and overlying shattered foam is located at the depth where the volume fraction of bubbles has become equal to the value corresponding to close packing without significant deformation, ∼0.65.Below this boundary, Equations (44)-( 50) apply.Above the boundary, the overlying material is assumed to have a bulk porosity of 30%, corresponding to a loose random packing of irregular fragments.
Parts (a)-(e) of Figure 15 show the resulting profiles of the variation of porosity with depth in a flow for magma water contents, n, of 100, 200, 300, 500, and 750 ppm.In each case the upper zone with an assumed constant porosity of 30% is the collapsed remains of what would have been the highly vesicular top of the flow if it had not been unstable.There is then an abrupt transition to the stable 65% porosity of the top of the gas-filled foam, and then a decrease in porosity down to the depth where H 2 O is just saturated in the magma.Below that depth no water exsolves and there is no vesicularity.We have ignored sulfur compounds in these calculations but assume that they behave in a similar way to the more abundant water.Note also that the vesicularity profiles of Figure 15 are dictated solely by the magma volatile content and not by the thickness of the flow.Thus, for n = 300 ppm (Figure 15(c)) the vesicular zone extends from a depth of 3.5 m to a depth of 4.9 m.If the flow containing this vesicular layer were 10 m thick, the layer would occupy the region extending from 35% to 49% of its depth.If the flow were 5 m thick, the vesicular zone would extend from 70% to 98% of its depth.If the flow were 4 m thick, the vesicular zone would extend from 87.5% of its depth to the base of the flow.Finally, the implication is that no stable laminar lava flow less than 3.5 m thick could be formed from magma releasing 300 ppm H 2 O as it was erupting.The Note.For a range of residual mass fractions, n, of gas in the magma rising into a lava lake in the vent, values are given for the final gas bubble diameters, f f , the pressures in the bubbles, P s , the maximum speed of the pyroclastic droplets produced as the bubbles burst, U p , and the maximum height reached by the droplets, H s .As lava travels away from the vent as a fluidized bed values are given for the now-reduced pyroclast ejection speed, V, the maximum height reached, H s , and the clast travel time, τ str .minimum stable flow thicknesses for n = 100, 200, 300, 500, and 750 ppm are 0.5, 1.7, 3.6, 8.6, and 16.9 m, respectively.We discuss the implications of this process of "auto-regolith" formation elsewhere (Head & Wilson 2020).Where stable flows of the kind illustrated in Figures 14 and 15 form, we stress that they do so late in an eruption (Phase 4; Figure 2) and hence form the proximal part of the resulting lava flow field.

Post-emplacement Modification of Lava Flows
The best-preserved of the Mare Imbrium flows terminate in the topographically lowest part of the mare-filled Imbrium impact basin (Schaber 1973).Although we have shown in Section 3 how mare flow length is limited by the internal rheological and thermal structure of the flow, it is still possible that the advance of the distal parts of these flows was modified by the topography-essentially, they ponded in the lowest part of the basin.In all cases, the strong possibility also exists that vesicular lava erupted late in the eruption, intruding the stationary distal flow deposits and inflating them (Phase 4; Figures 2(a), (e)), both by simple injection of additional lava (Figure 12) and by additional vesiculation of the lava during cooling-induced second boiling.Just such a sequence of emplacement has been inferred for flood basalt flows on Earth (Self et al. 1996(Self et al. , 1997;;Thordarson & Self 1998), where inflation by at least a factor of 3 and as much as a factor of 10 has been suggested.We now explore how modification processes may have acted in various parts of an emplaced flow deposit.

Proximal Modification
Where activity has already continued for days to weeks (Phase 4; Figures 2(a), (e)), relaxation of the dike induced by tectonic stresses is accompanied by cooling of the magma against the dike walls and a general decrease in the volume flux of partly vesicular magma emerging from the vent into a lava lake and proximal flows.Lava in the outer part of the pond around the vent may have developed a significant cooled crust due to a decrease in the radius of the fire fountain as more of the volatiles are included in Strombolian explosions and less into steady magma fragmentation.Deformation of this crust as a result of the effusion rate decrease can cause fractures to form, through which hot vesicular magma may emerge.These processes have been invoked by Qiao et al. (2017Qiao et al. ( , 2018bQiao et al. ( , 2020aQiao et al. ( , 2020b) ) to explain low dome-like structures, seen on the floor of the small shield volcano crater lake at Ina, in the volcanic depression Sosigenes, and in the small shield volcano Cauchy 5 summit pit (Qiao et al. 2018a(Qiao et al. , 2020a)).These are one subclass of the features described as irregular mare patches (IMPs) by Braden et al. (2014).
A much more dramatic form of modification of laminar proximal flows emplaced late in the eruption is likely to be common.Ongoing crystallization due to cooling during and after emplacement will cause supersaturation of the volatiles remaining in the lava and induce second boiling.The low pressures dictated by the absence of atmosphere and low acceleration due to gravity on the Moon mean that this process will operate much more effectively than on Earth.Volatile exsolution does not all take place at the same time.It is driven by the penetration of waves of cooling propagating from the top and base of the flow and so the central part of the flow, consisting of the existing foam layer and the unvesiculated lava beneath it (right-hand part of Figure 14) is the last to be affected.
Consider Figure 15(b), where the emplaced proximal flow contained 200 ppm water.In the region occupied by foam, the vesicularity decreased from 65% to zero over the depth range 1.75 to 2.18 m.The vesicularity was zero at all greater depths, with the lava retaining 183 ppm water at 1.75 m depth increasing to 200 ppm retained at all depths greater than 2.18 m.Forced exsolution of all the residual water would increase the vesicularity of the lava to 95% at 1.75 m depth and 93% at 2.8 m depth, inflating the flow in the process.The vesicularity would have reached 90% at a depth of ∼27 m, the lava at this level having originally been located 3.7 m below the surface of the flow.The implication is that if the flow had originally been 3.7 m thick, its thickness would have potentially grown to ∼27 m, a more than sevenfold increase in thickness, if the foam of which it now almost entirely consisted had remained stable.It is, however, very unlikely that it would have remained stable during this expansion process.Thermal stresses during cooling or, eventually, shock waves from the impacts of micrometeorites are likely to have caused fracturing of the crust, allowing disruption of the expanded foam.Escape of the liberated gas through fractures in the 1.75 m thick crust would have allowed the crust to subside, while the expanding gas carried the fragments of collapsed bubble walls upward to eventually settle back over the subsided slabs of crust.Proximal flow disintegration driven by this process may be the origin of a second type of IMP such as those seen at the small shield volcano Cauchy 5 (Qiao et al. 2018a) and at numerous other locations in mare lava flow fields (Wilson & Head 2017b;Qiao et al. 2020b), and may also contribute to the formation of auto-regolith (Head & Wilson 2020).

Medial and Distal Modification
Section 3.3 showed that, during their actively advancing stage, the medial and distal parts of large-scale lunar flows will consist of lava that has lost most of its volatiles, both in the lava fountain and during the turbulent flow phase through which it has passed.The lava will therefore be largely vesicle free.In addition, the plug formed as the lava motion changes from turbulent to laminar will occupy the upper ∼75% of the flow (Figure 9(h)).As a result, if distal intrusion and inflation occur after the flow front ceases to advance, intrusion will likely occur just below the base of the plug.The intruding lava may possess a vesicular core, similar to that emplaced close to the vent (Figure 15), but the vesicular layer thickness will have been significantly reduced during the intrusion process because the density, and hence weight, of the plug lava will have been greater than the fragmental layer density used to calculate the profiles in Figure 15.Furthermore, as was the case for the proximal lava, the late-stage lava may not only have already contained a vesicular layer, but may also have consisted, in part, of lava that had not had the opportunity to release all its volatiles.In this case, as shown in Section 4.1 for proximal flows, cooling-driven second boiling has the potential to inflate all parts of the flow field that have been intruded.The extent of this process can again be illustrated using Equations ( 44)-( 50), where now the same spreadsheet used for the earlier calculations is extended by replacing the shattered upper layer of Figures 14 and 15 with the largely vesicle-free plug and allowing the intruded magma to exsolve volatiles as it crystallizes.The same issues arise as for proximal flows: It is easy to release enough gas through second boiling to dramatically inflate a flow until its interior foam layer becomes unstable and collapses.
Figure 16 shows the results.The solid lines represent five flows with initial thicknesses of 2, 4, 6, 8, and 10 m.Each flow is intruded by lava forming a layer whose thickness is 20% of the flow thickness.This corresponds to the intrusion taking place in a flow whose front has ceased to advance after 80% of the total lava volume has been erupted.The top of the intruded layer is located at a depth below the flow surface that is 75% of the flow thickness, placing it at the base of the plug.The solid lines show the final thickness of each of the five flows after the complete exsolution of from zero to 60 ppm of dissolved water.Note that even with zero water release the flow inflates by 20% solely due to the intrusion.When water release does occur, the final state of the lava is controlled by the initial flow thickness and the amount of gas exsolved.The broken lines in Figure 16 define the combinations of final thickness and exsolved water at which the bubbles of gas would occupy 65%, 75%, and 85% of the volume of the intrusion, and therefore represent likely stability criteria.Any final flow thickness above a given dashed line would avoid collapse and represents a stable inflated flow (at least until regolith development begins; Head & Wilson 2020).For the conservative stability limit of 65% gas bubble content used for the proximal flows, an initially 2 m thick flow could inflate to 3.16 m and remain stable.If the stability criterion were relaxed to 75%, the same flow could inflate to just over 3.6 m.With the 65% stability limit, an initially 10 m thick flow could inflate to ∼15.8 m and remain stable.If the stability criterion were relaxed to 75%, the same flow could inflate to just over 18 m, and if the stability limit were relaxed to 85% it could inflate to ∼26 m before collapse.In summary, this process can readily lead to stable distal flow deposits inflated by a factor close to 2, with ∼50% of the flow consisting of very vesicular lava.
Figure 17 shows an explicit example of the successive developments (left to right in the figure) during the inflation process for a large-volume mare lava flow.At the end of its laminar flow stage, it comes to rest with a thickness of 8 m.The hot lava core is vesicle free, having lost all its gas in the fire fountain at the vent, and has a crystal content of 30%.The solidified upper and lower thermal boundary layers are ∼1.3 m thick.This flow is now intruded by lava from the still-active vent such that it inflates by 1.6 m.Assuming that the inflation takes place throughout the flow, this implies that 20% of the total lava volume is erupted from the vent after the flow front stops advancing.The intruded lava is assumed to have a residual water content of 100 ppm so that, after mixing with the hot lava already present, the average water content of the liquid is 30 ppm.During the following 30 days, cooling fronts propagate more than half of the way to the center of the flow, causing additional crystallization and concentrating the water into the remaining liquid lava so that supersaturation occurs.The resulting gas bubbles inflate the lava so that the flow thickness is now 13.8 m.The process continues until, by 55 days after emplacement, the entire depth of the flow is partly crystallized, and the thickness has increased to 17.7 m, more than twice the original value.During this period, various processes of gas bubble segregation or incipient convection may have occurred.The upper crust of the flow will have experienced stresses due to cooling and to initial flow thickness variations caused by preexisting topgraphy, which will have led to variations in the amount of inflation from place to place.These stresses may crack the crust, and the final part of Figure 17 shows two possible consequences: the expulsion of very viscous foam to form ring-moat dome structures (RMDSs; Zhang et al. 2017), or the partial disintegration of the upper part of the flow, a possible explanation of IMPs (Braden et al. 2014).
Estimates in the literature of the volumes of individual lunar lava flows, used in flow emplacement modeling, and of the total volumes of lava present in the various mare basins, which have implications for lunar thermal history, have generally assumed that the observed or inferred volume is the dense rock equivalent volume.If inflation and second boiling are common in lunar lava flows, these volumes may be overestimated by at least as much as a factor of 2. Additionally, indirect thickness and volume estimates based on, for example, radar sounding, derive flow thicknesses by assuming a density of the substrate, therefore such analyses may themselves be subject to systematic errors.

Discussion
We have applied the basic principles of magma generation, ascent, and eruption derived from analysis of terrestrial basaltic eruptions and developed a theoretical treatment of the dynamics of lunar mare basalt volatile release and lava flow emplacement and evolution.Unlike mafic magmas on Earth (Figure 1), lunar basaltic magmas derive from deep mantle magmatic sources, encounter global low-density crust, are erupted into a low-gravity environment characterized by lack of an atmosphere, and are one to three orders of magnitude less viscous than terrestrial basalts (Rai et al. 2019, and references therein).Exploration of the consequences of the different environments of terrestrial and lunar eruptions (Figure 1) forms the basis for a new paradigm, one that often seems counterintuitive when viewed from our traditional Earth-based experience.These factors combine (Figure 2) to make lunar lava flows significantly different from their terrestrial counterparts in their (i) total volume, (ii) magma rise speed and volume flux, (iii) stages of evolution, (iv) vesicularity, and (v) resulting vertical and horizontal structure.We identified four phases in lunar eruptions and lava flow emplacement: Phase 1, initial gas venting phase.Initial penetration of a deep-sourced dike to the surface and dike-tip gas venting.
Phase 2, Hawaiian eruptive phase.High-flux (10 6 -10 5 m 3 s −1 ) explosive Hawaiian eruptive phase, with ∼70% of the dike volume emerging as a largely degassed lava flow from a pyroclast-fed lava lake and flowing at ∼5 m s −1 for up to about 300 km, initially as a turbulent flow and then as a laminar flow.
Phase 4, Strombolian vesicular flow phase.Finally, dike closure forces extrusion of the remaining ∼25% of the dike volume as a low-flux (∼10 4 m 3 s −1 ), low-speed (∼0.1 m s −1 ) laminar and more vesicular flow.During this final phase, flow front advance eventually ceases due to a combination of rheological and thermal effects (Figure 8(e)), and the stillmolten core of the cooling flow can be intruded by volatilecontaining magma from the source region (Figure 14).Flow inflation can occur, and as the flow continues to cool and solidify, second boiling can produce very vesicular magmatic foam layers that can potentially extrude through cracks in the overlying flow (Figure 17).
The documentation of these four basic lunar lava flow eruption phases provides the opportunity to explore the consequences for volatile release patterns, flow morphology, vertical and lateral flow structure, the consequent development of regolith, and to propose observational tests to further develop and refine models of lunar basaltic volcanic eruptions and lava flow emplacement.

Volatile Release Patterns and Surface and Near-surface Structure of Proximal Flows
Phase 1 of a lunar eruption (Section 3.2) involves the escape of a quantity of gas that has accumulated immediately below the upper tip of a rising dike as the tip reaches the surface.This phase is extremely explosive as free gas in the tip of the rising dike is erupted when the dike tip reaches the surface.Pyroclasts are distributed over a large part of the lunar surface, forming a thin, undetectable deposit.However, this phase lasts only a few minutes and involves only ∼0.25% of the total magma mass erupted and has no bearing on the nature of lava flows and pyroclasts produced subsequently (Wilson & Head 2017a).
Eruption Phase 2 (Section 3.3) lasts tens of hours and is also explosive, forming a lava fountain.The shape and dynamics of the fountain (Morgan et al. 2021) have much in common with the largest volcanic plumes seen on Io, but on a much smaller spatial scale, because the intrinsic volatile contents of lunar magmas are much less than those of Io magmas, which entrain sulfur compounds as they rise through Io's crust (Leone & Wilson 2001).The heights of lunar lava fountains range from 875 m to 8.75 km for total volatile contents in the range 300 to 3000 ppm, with corresponding widths ranging from 3.5 to 35 km (Morgan et al. 2021).The small sizes of lunar pyroclasts coupled with the large volume fluxes expected during Phase 2 ensure that lava fountains are generally optically dense, so that most of the clasts cannot cool during flight.Decompression of the expanding volcanic volatiles to the lunar vacuum ensures that pyroclasts lose most of their volatiles in flight, so that the falling degassed droplets accumulate into a very hot lava lake around the vent.The lava flows leaving this lake have extremely low vesicularities.The progressively decreasing volume flux during the latter part of Phase 2 of an eruption causes an increasing fraction of the outer shell of the lava fountain to become translucent, so that the pyroclasts in that part of the fountain can radiate heat.Clasts landing at any given distance from the vent can reach that range via a high-angle or a low-angle trajectory, and so cooling of the outer shell of the fountain causes a general decrease in the temperature of the pond created by the falling clasts, with an additional bias toward the clasts falling at the greatest range being the coolest.This can promote the accumulation of clasts into a spatter or cinder cone deposit at the margin of the lava pond.Based on likely feeding dike geometries (Wilson & Head 2017a), ∼70% of the magma volatiles, together with 70% of the dike magma (Figure 2(a)) are likely erupted during Phase 2 of the eruption.
The decreasing volume eruption rate that marks Phase 3 (Section 3.4) is the result of more of the rising dike feeding the eruption being located in the crust, where magma is negatively buoyant, than in the mantle, where it is positively buoyant.Reduced magma rise speed allows time for the coalescence of exsolved gas bubbles, especially the CO bubbles nucleated in the deepest parts of the dike.Coalescence produces a progressively greater proportion of unusually large bubbles in the magma reaching the surface.By the time the largest of these are big enough to occupy a large fraction of the width of the dike, their intermittent emergence though the vent causes the lava fountain to become increasingly unsteady and pulsatory.Furthermore, the largest gas bubbles do not interact efficiently with the pyroclasts being formed by the fragmentation of the rest of the magma, and so the gas mass fraction providing the energy driving the explosive process decreases.As a result, the lava fountain height and width begin to decrease.
Eventually the lower tip of the dike feeding the eruption ceases to rise as the dike as a whole reaches buoyancy equilibrium with the density structure of the crust and upper mantle.This initiates Phase 4 (Section 3.4).The only process driving magma out of the dike is now the slow reduction of the horizontal width of the dike as the stresses in the lithosphere caused by the presence of the dike relax.The very low magma rise speed allows more time for gas bubble coalescence, and eventually the largest gas bubbles occupy almost all the width of the dike.They are then classed as slugs, and their emergence though the vent produces intermittent Strombolian explosions.These take place simultaneously with the explosive ejection of the magma through which the slugs emerge.Even greater scavenging of bubbles of water and sulfur compounds by the slugs causes the effective gas mass fraction driving the fountain to become even smaller.This causes the fountain height and width to decrease further.Equations ( 7) and (8) show that a reduction in fountain radius, R, causes a decrease in the extent, X, of the outer shell through which heat escapes, and so if the decrease in R is greater than the decrease in the volume flux, F, of magma being erupted, it is possible to have a transient situation where the lava fountain gets smaller but hotter.Nevertheless, there must come a time late in Phase 4 when the residual magma in the dike is very volatile depleted, and the lava erupted between the Strombolian explosions contains only a few hundred parts per million of water and sulfur compounds.Fragmentation of this magma produces pyroclastic droplets, forming lava fountains that are at most several tens to a few hundred meters high immediately above the vent (Section 3.4; Table 4).Interference between ascending and descending clasts decreases the fountain height further with increasing distance from the vent, and the upper part of the lava flowing away resembles a fluidized bed.As the remaining gas escapes, the bed becomes compacted and the nature of the last lava flowing away from the vent lava pond is predicted to be a layer of loose fragments, overlying a layer of welded fragments, in turn overlying a shearing foam lava flow with negligible vesicularity at its base (Figure 14).The remaining ∼30% of the magma and its associated volatiles are erupted during this last phase of the eruption.

Implications for Regolith Formation and Evolution: Initial
Development of "Auto-regolith" Most mare regolith development models begin with the assumption that solid basaltic bedrock is the initial surface undergoing impact bombardment, energy partitioning, mechanical fragmentation, and soil (regolith) development (e.g., McKay et al. 1991;Spudis & Pieters 1991, their Figures 10-13;Wilcox et al. 2005, their Figure 12).It is clear from our analysis that the final stages of lava flow emplacement may result in a surface that is very different from solidified basaltic bedrock, and that the flow interior may also differ substantially from solid bedrock (Figure 15).For example, proximal regions of typical lunar lava flows like the Phase 2 and 3 flows in Mare Imbrium (Schaber 1973) may be characterized by fresh surfaces dominated by several-meter-thick deposits of glass shards from exploding vesicles, overlying highly vesicular lowdensity layers (Figure 15), before transitioning down into highdensity basaltic bedrock.This configuration will have a significant influence on energy partitioning in initial impact events, with relatively more energy going into glass shard redistribution and bubble/vesicle wall crushing, rather than into brittle deformation of solid basalt (Head & Wilson 2020).In addition, the internal structure of the flow will have an influence on regolith development: Penetration of impacts into internal vesicle/foam layers will complicate simple models of transient cavity development, energy partitioning, and final fresh crater form.This configuration will have very significant effects on the regolith thickness estimates based on the assumption of regolith initiation in a solid basaltic substrate surface (Head & Wilson 2020).In the cases where a fragmental layer is initially produced during lava flow emplacement, the regolith thicknesses generated by impact processes will be significantly overestimated.We term this process auto-regolith formation due to the fact that the nature of the eruption, degassing and solidification of the lava itself produces an initial surface fragmental soil layer that could be of substantial thickness, making up a large part of the regolith thickness observed today, and altering the perception of the thickness of impact-induced regolith (for a summary of the implications of auto-regolith, see Head & Wilson 2020).Furthermore, for areas characterized by RMDS (Zhang et al. 2020(Zhang et al. , 2021)), interpreted to be magmatic foams, impact event energy partitioning should be dominated by crushing and penetration (the aerogel effect), as opposed to brittle deformation and lateral ejection, resulting in smaller, deeper craters and less ejecta.Therefore, this effect may also alter the superposed impact crater size-frequency distribution, yielding anomalously younger ages and potentially invalidating the application of current diffusional landform degradation models.

Vertical Structure of Lava Flows
Relatively thin (<∼10 m) terrestrial lava flows typically display a vertical structure composed of three components (Aubele et al. 1988): (i) an upper vesicular zone (about one half of the vertical section of the flow; vesicles decrease in diameter and increase in number toward the top), (ii) a middle nonvesicular or dense zone, and (iii) a lower vesicular zone (about 30-40 cm thick regardless of flow thickness; vesicles increase in diameter and decrease in number toward the top of the zone; Figure 1(f)).Aubele et al. (1988) interpret these zones to be formed due to volatile exsolution, growth and rise of bubbles during cooling, and flow solidification subsequent to flow emplacement.This tripartite flow structure, as well as vesicle density and size, are interpreted to be related to the advancing upper and lower lava solidification fronts.Lunar lava flows will undergo similar volatile release, vesicle formation and cooling processes, but the details will be different due to the unusual lunar environment (Figures 14,15,and § 17).The low lunar gravity and lunar vacuum will cause the upper vesicular layer to be much thicker and more complex, with bubble sizes being larger than on Earth, and the lower vesicular layer should have even larger vesicles on average.The lack of atmosphere and hence of convective cooling of flow surfaces will influence boundary layer growth rates and thicknesses.Furthermore, the flow interior can be made more complex by the process and aftermath of flow inflation (Figure 17), the injection of new lava into the cores of just-emplaced flows, thought to have been typical of many terrestrial flood basalts.If these inflated flows undergo second boiling upon cooling and solidification, additional vesicle layers can be formed.The observed terrestrial and predicted lunar flow cross sections (compare Figures 1(f) and 17) illustrate these different consequences and could help interpret the exposures of vertical sections of lava flows, as seen in lunar crater and rille walls that are candidates for future lunar access and exploration.

The Role of Second Boiling in Lava Flow Emplacement
We have shown that, during Phase 4, late-stage lava inflating distal parts of flows would have contained a vesicular core; both the distal and proximal parts of the flow field could have contained late-stage vesicular lava (Figure 17).In addition to containing a vesicular layer, this lava may have consisted, in part, of lava that had not yet had the opportunity to release all its gas.After the eruption ended, if that volatile-containing magma flow component underwent second boiling, i.e., was forced to release volatiles as it continued to crystallize, then the flow field will have undergone inflation, increasing in thickness by a factor of up to ∼3 (Figure 17).Gas exsolution does not all take place at the same time, beginning both above the base of the flow and beneath the top of the flow as the two waves of cooling extend into the flow, and ever more new bubbles nucleate as the supersaturation increases, until all the dissolved volatiles have been released.The last place to undergo this vesiculation is close to the middle of the flow.If fractures occur in the cooled crust of this flow, hot vesicular lava foam from the interior can extrude onto the surface of the flow, creating low mounds and surrounding depressions.This is a suggested mechanism (Wilson et al. 2019) for the formation (Figure 17) of distinctive low mounds on the surfaces of mare flows classified as RMDSs by Zhang et al. (2017).

Importance of Lunar Lava Flow Predictions for Future Exploration
The predictions outlined in this paper provide an important framework for the future documentation and characterization of lunar lava flows by remote-sensing methods and human-robotic exploration.Although lunar mare lava flows are currently covered with regolith to several meters depth, smaller regolith thicknesses will have accumulated between flows formed earlier in lunar history when the time intervals between eruptive events were shorter.Thus, many of the original characteristics of the flow phases will have been preserved.These properties can be explored with remote radiometer, spectral reflectance and radar data, and in situ analyses, including outcrop observations at skylights and at crater and rille walls, and with ground-penetrating radar (e.g., Fa et al. 2015;Lai et al. 2016).

Importance for Interpretation of Past Exploration Results
and Planning for Future Experiment Design During the Apollo Lunar Exploration Program, a variety of surface experiments were deployed (seismic, heat flow, surface electrical properties, gravity) that required models of the regolith and bedrock structure for data interpretation.Similarly, orbital experiments (e.g., Apollo lunar radar sounder, Kaguya radar, Lunar Reconnaissance Orbiter Diviner, etc.) also utilize models of near-surface physical properties for interpretations.Finally, surface-rover exploration experiments (e.g., Chang'e ground-penetrating radar) rely on assumptions about regolith and bedrock density and porosity.Documentation of the predicted range of properties of lunar lava flows that have passed through the four phases we have proposed (e.g., Head & Wilson 2020) may help in the interpretation of data from these experiments and can be factored into the design and interpretation of future mission experiments.

Conclusions
Lunar mare basalt eruptions followed a predictable pattern involving four phases (Figure 1).After a very-short-lived ventopening and gas venting Phase 1, a Hawaiian fire-fountain Phase 2 initiated flow formation.The measured low viscosities and predicted high initial eruption rates of lunar mare basalts imply that the motion of lava leaving the vent was commonly turbulent, despite the consequent efficient radiative cooling of the bulk of the lava leading to the rapid development of non-Newtonian rheology (Figures 8 and 9).
For lavas that went through a turbulent phase on eruption, the turbulence persisted for a distance from the vent that was typically 20%-25% of the final length of the flow (Figure 10).Further, the absolute distance from the vent at which the turbulent-to-laminar transition took place depended almost entirely on the volume flux per unit width of the flow (Figure 10(a)).This result has great significance for eruptions that continued long enough for turbulent heat transfer at the base of the flow to have caused thermal erosion of the substrate, leading to sinuous rille channel formation.Measuring the length and mean width of the channel immediately provides an estimate of the volume eruption rate.
For flows with an initial turbulent phase, the initial ∼40% (Figure 11) of cooling-induced crystal growth occurred throughout the flow at rates of ∼10-20 K hr −1 (Figure 13).The remaining 60% of crystallization occurred during the subsequent laminar phase at rates from one to several orders of magnitude slower, the rate depending on a sample's final depth in the flow deposit.The crystals in any given sample of lava from such a flow should therefore show distinctive radial variations of cooling rate.
After a transition (Phase 3) to low-effusion-rate Strombolian activity (Phase 4), syn-eruptive flow inflation (Figures 14 and  15) due solely to the eruption continuing after the flow front stopped advancing was commonly by a factor of 1.5-2.5 (Figure 12).Further inflation by an additional factor of up to 2 (Figure 16) was possible after flow front advance ceased if second boiling drove residual volatiles out of solution as the lava cooled to form a very vesicular foam.Escape of foam through fractures in the crust may have formed RMDSs.In some cases, cooling stresses in the resulting vesicular lavas (Figure 17) may have led to foam collapse and subsequent partial re-compaction, forming IMPs.
The above predictions form a framework for the interpretation of existing data on mare lava flows and the planning of future lunar mission experiments.

Figure 1 .
Figure 1.Comparison of magmatic-volcanic-tectonic-surface environments on the Earth and Moon.Differences are shown in global geodynamic environment; crustal type, thickness and continuity; depth of magma origin; ascent and eruption (effusive and explosive); surface environment; flow cooling behavior; type of flow; and flow structure.(a).The range of magmatic-volcanic environments on the Earth.(b).Characteristics of explosive and effusive environments on the Earth.(c).The range of magmatic-volcanic environments on the Moon.(d).Comparison of magmatic volatile release patterns on the Earth and the Moon.(e).Terrestrial lava flow vent structure and flow elements.(f) Generalized structure of terrestrial flood-basalt-scale lava flow synthesized from Aubele et al. (1988), Self et al. (1996, 1997), and Thordarson & Self (1998).

Figure 3 .
Figure 3. Variation of mean magma rise speed, U d , at great depth with horizontal outcrop length, L, of erupting fissure.Values are given for a series of dense rock equivalent volume fluxes, F, of magma being erupted (lines labeled in cubic meters per second).The shaded area marks the transition from Hawaiian to Strombolian explosive activity as the erupted volume flux decreases during an eruption.The upper edge of the shaded region corresponds to magmas rich in CO (∼1400 ppm) and the lower edge to CO-poor magmas (<100 ppm CO).
-source vent forming a circularly symmetric fountain (Figure 4(a)).For a fissure vent erupting actively for a distance L along its strike (Figure 4(b))

Figure 4 .
Figure 4. Relationships controlled by the relative values of the range, R, of pyroclasts and the length, L, of the fissure that erupts them.(a) If R is much greater than L, a near-circular pond forms around the vent and the number density of pyroclasts in the fire fountain decreases as the square of the distance from the vent.(b) If L is much greater than R, an elongate pond forms and clast number density decreases linearly with distance, except at the ends of the fissure.
Pinkerton & Stevenson (1992) for larger crystal fractions.Next, we evaluate the yield strength of the lava.Functions based on experimental data approximating the variations of yield strengths of fluids with their solids contents are given byGay et al. (1969),Ryerson et al. (1988),Dragoni (1989), Zhou et al. (1995), Saar et al. (2001), Mueller et al. (2010), Pantet et al. (2010), and Ishibashi & Sato 3.3.3.Implications for Dynamics and Morphology of Lava Leaving the Vent PondLava rheology depends on the crystal content f c , which is a function of the lava temperature T e via Equation (12a), in turn dependent on the hot pyroclast fraction f hot via Equation (11), and ultimately determined by the pyroclast range R and the erupted volume flux F, via Equations (8) and (9).Using these equations, we show in Figure7the relationship between f c , R, and F.
Figure 7(a) is for an eruption creating an elongate pond from a fissure vent 10 km long, a typical value for mare lava fissure vents (Head & Wilson 2017).
Figure 7(b) is for an eruption from an extremely short fissure, effectively a pointsource vent, so that the lava pond is approximately circular.

Figure 6 .
Figure 6.Variation of yield strength Y (left ordinate) and viscosity correction factor f v (right ordinate) with lava crystal content f c for any mafic lava.Also shown, using the right ordinate scale, is the plastic viscosity η p of lava with the composition of the Apollo 12 basalt sample 12002 (see Williams et al. 2000).

Figure 5 .
Figure5.Variation of the percentage of a lava lake around a fissure vent that forms from pyroclasts at magmatic temperature.The fraction of pyroclasts remaining hot in a dense lava fountain is a function of the maximum range of the pyroclasts, which defines the maximum extent of the lake from the vent, and the total dense rock equivalent magma volume flux being erupted from the vent (curves labeled in cubic meters per second).Upper abscissa shows the released magma volatile contents (dominantly CO) that correspond to the ranges on the lower abscissa.(a) Fissure vent with outcrop length of 10 km, (b) fissure length of 3 km.

Figure 7 .
Figure 7. Variation of the volume% of crystals f c in lava leaving a vent pond (left ordinate) with erupted magma volume flux F and maximum pyroclast range R (labeled lines) defining the extent of the pond.(a) Elongated pond from a 10 km long fissure, (b) near-circular pond from a very short fissure.In each case the right ordinate shows the rheological parameters corresponding to f c from Figure 6.
) and (c) show plan views and Figures 8(b) and (d) show cross sections of non-Newtonian flows of lava with Bingham plastic rheology.Figures 8(a) and (b) correspond to laminar flow and Figures 8(c) and (d) correspond to turbulent flow.In Figures 8(a) and (c), stationary levees are indicated at the edges of the flows, and in all four figures the lava velocity profiles are shown using broken curves relative to an axis.In laminar flow, an unsheared layer called a plug is carried along by the shearing lava beneath it (Figure 8(b)), which has a parabolic velocity profile (Figures 8(a) and (b)).In turbulent flow (Figures 8(c) and (d)) the lava is fully mixed, producing a less strongly curved flow front (compare Figures 8(a) and (c)) and so-called "top-hat" velocity profiles (Figures 8(c) and (d)).

Figure 8 .
Figure 8. Characteristics of non-Newtonian (specifically Bingham) lava flows.(a) and (c): plan views, (b) and (d): vertical axial cross sections.(a) and (b): turbulent flow, (c) and (d): laminar flow.Velocity profiles are superposed as short dashed lines.(e): representation of the internal structure of the flow as the change from turbulent to laminar occurs.Thermal boundary layers grow only after turbulence ceases.Lava motion is toward the top in (a) and (c), and toward the right in (b), (d), and (e).
(27) reduces to the Jeffreys equation for laminar Newtonian fluids, U l = [(ρ g D l 2 sin α)/ (3 η)].In this case, we can use Equation (21) to eliminate U l from Equation (27) and find When the top of the basal thermal boundary layer rises to meet the bottom of the plug, i.e., when h + z = D f , no further shearing can occur anywhere within the flow, and the flow front must stop (Figure 8(e)).If the vent is still erupting, then either a breakout occurs somewhere (see Head & Wilson 2017, their Figure 12(b)) (a)), the exponential flux decay constant τ d (Figure 10(b)), and the surface slope α (Figure 10(c)).In each case the other two variables are held constant at intermediate values.The maximum flow length increases with increasing values of all three variables F/L, τ d , and α.It is clear from Figure 10, derived from Table

Figure 9 .
Figure 9. Variation with distance from the vent of 10 parameters of the initially turbulent lava flow listed as model #1 in Table 2. (a) Mean flow speed, (b) flow thickness, (c) absolute temperature of the entire flow thickness while turbulent and of the flow core between growing thermal boundary layers while laminar, (d) volume fraction of the flow consisting of solids, (e) lava yield strength, (f) lava plastic viscosity, (g) thickness of the lower thermal boundary layer, (h) thickness of the parts of the flow that cannot be sheared due to cooling or the presence of a yield strength, (i) thickness of the part of the flow that consists of lava that can be sheared, (j) rheological parameters Reynolds number, Hedström number, and critical Reynolds number.

Figure 10 .
Figure 10.Variation of the maximum length of a flow (upper curve) and the distance from the vent at which turbulence dies out (lower curve).Both are shown as a function of (a) the erupted volume flux divided by the width of the flow F/L, (b) the exponential decay time of the eruption rate τ d , and (c) the slope of the surface α on which the flow travels.

Figure 11 .
Figure 11.Total solids volume fraction in a flow at the point when the transition from turbulent to laminar flow occurs.Values are shown as a function of (a) the erupted volume flux divided by the width of the flow F/L, (b) the exponential decay time of the eruption rate τ d , (c) the slope of the surface α on which the flow travels.

Figure 12 .
Figure 12.Factor I by which the thickness of a flow increases after the flow front ceases to advance and the body of the flow inflates due to continuing activity at the vent, as a function of the erupted volume flux divided by the width of the flow F/L.The continuous curve is for an exponential flux decay time τ d of 1 day and a ground slope α of 0.002.The two short segments show the effect of increasing τ d to 2 and then 3 days.The just-visible vertical line extending above and below the α = 0.002 curve at F/L = 30 m 2 s −1 demonstrates the very small influence of α varying between 0.001 and 0.004.

Figure 13 .
Figure 13.Average cooling rate C t of lava flowing away from a vent in a turbulent manner as a function of the erupted volume flux divided by the width of the flow F/L.The vertical scatter of data points at F/L = 30 m 2 s −1 represents the range for the data in Table2, i.e., α between 0.0005 and 0.004, and τ d between 0.5 and 3 days.

Figure 14 .
Figure 14.Development of flow structure in the immediate vicinity of the vent as volatile-containing magma erupts into the lunar surface vacuum environment during Phases 3 and 4 as the transition from steady Hawaiian to fully Strombolian activity occurs.Gas slugs rise through the conduit and burst through the lava surface to cause Strombolian explosions, propelling clots of viscous magma through a fluidized bed generated by the explosive decompression and disruption of very vesicular lava forming a foam layer at the top of the magma column in the vent.The final products as lava moves away to the right are a surface layer of loose glass shards overlying partly to fully welded pyroclastic fragments, in turn overlying a shearing foam layer, and finally a basal layer of nonvesicular lava that still contains volatiles that can later be driven out of solution by second boiling as lava crystallizes.

Figure 15 .
Figure 15.Variation of vesicularity with depth in lava flows formed during the fully Strombolian Phase 4 of an eruption.Profiles are shown for magma water contents n of (a) 100, (b) 200, (c) 300, (d) 500, and (e) 750 ppm.The vertical thickness of the flow below the depth where the vesicularity goes to zero does not influence the patterns of vesicularity shown.Frame (f) gives a visual impression of the basal unvesiculated lava, overlain by the vesicular layer, in turn overlain by the accumulated fragments from the explosive decompression of the underlying foam.The image of the fragments is a USGS photograph of lapilli from the Pu'u Pua'i eruption, Kilauea, Hawai'i.

Figure 16 .
Figure 16.Final thicknesses of inflated flows with initial thicknesses of 2, 4, 6, 8, and 10 m (labeled solid lines) as a function of the mass fraction of H 2 O exsolved to form a foam during crystallization-induced second boiling of intruded lava.The intruded lava represents 20% of the initial flow thickness in each case.Broken lines define combinations of final thickness and exsolved H 2 O at which the gas bubbles occupy 65%, 75%, and 85% of the volume of the intrusion, values representing possible foam stability criteria.

Figure 17 .
Figure 17.Five stages in the evolution of the internal structure of a distal flow during the final stages of flow emplacement, at final rest, and following further flow cooling and second boiling for the cases where late-stage lava flows have ∼100 ppm volatile content.(a) Internal structure of distal flow during late stage of flow emplacement and at final rest.(b) Flow inflation phase.(c) Flow cooling, crystallization, and second boiling phase (after 30 days).(d) Flow cooling, crystallization, and second boiling phase (after 55 days).(e) Possible consequences of late stage of second boiling: vesicular lava (foam).The cold basalt crust may fracture due to accumulating stresses allowing foam to escape through cracks to form a ring-moat dome structure (RMDS).Alternatively, the fractured crust may subside through the foam allowing explosive decompression of the top of the foam to form a more compacted fragmental layer, analogous to part (f) of Figure 15.Net depressions formed in this way in otherwise intact parts of a flow may be the origin of irregular mare patches (IMPs).This figure is based on Figure 5 in Wilson et al. 2019) but is calculated for a flow thickness similar to that shown in Figure 9 of this paper.It has been modified in part (d) to take account of the collapse of extremely vesicular lava allowing gas escape and extended in part (e) to include IMP formation.

Table 1
Consequence of Turbulent Lava Flows Entraining the Regolith Over Which They Flow Lava erupted at the liquidus temperature 1713 K entraining a volume fraction f r of loose clasts as it advances acquires a bulk temperature T e and a total volume fraction of solids f s as the clasts partially melt and the lava partially crystallizes.
(Skelland 1967)ained, we can find U t from Equation (21), Re from Equation (20), and f f from Equation (19).It is a requirement for turbulent flow that Re must be greater than some critical value, Re crit .For Newtonian fluids Re crit is commonly taken as a constant with value ∼2200, but for Bingham plastics it is not a constant but instead is a function of a second dimensionless number, the Hedström number for the flow, He, defined by(Skelland 1967) The lava temperature is 1592 K, Λ is 11.032 4, τ m /τ f is 72.63, and melting would not start until 26.2 days after the eruption begins, and only then if the discharge rate remained constant.At an intermediate distance of 50 km from the vent, first reached by lava after 2.93 hr with a temperature of 1660 K, Λ is 7.867, τ m /τ f is 38.62, and melting starts after 4.7 days if the eruption rate does not change.

Table 2
Parameters from 17 Model Eruptions Model

Table 3
Initial Lava Cooling RatesNote.The cooling rate, C e , varies strongly with depth below the flow surface.

Table 4
Conditions during Strombolian Activity Maximum ejection height of pyroclasts in Strombolian phase I Factor by which flow thickness increases due to inflation L Horizontal length of dike, equal to length of fissure erupting magma from dike L h Latent heat of fusion of magma, 6.06 × 10 5 J kg −1 Ratio (T sol − T a )/(T − T a ) T sol Solidus temperature of magma U d Mean magma rise speed in dike U f General value of lava flow speed U l Mean speed of laminar lava flow U p Pyroclast eruption speed U t Mean speed of turbulent lava flow U trans Mean lava flow speed where turbulent-to-laminar transition occurs V Volume of lava erupted in time t V tot Total volume of lava erupted V inf Volume of lava erupted when flow front stops and inflation starts W Mean width of dike X Thickness of shell in lava fountain from which clasts can radiate heat Y Lava yield strength Y trans Lava yield strength when turbulent-to-laminar transition occurs Z trans Distance from vent to location of turbulent-to-laminar transition Z max Total length of lava flow c c Mean specific heat of lava at subsolidus temperatures, 1000 J kg −1 K −1 c h Mean specific heat of lava in the melting range, 1570 J kg −1 K −1 d Mean diameter of pyroclasts, 300 μm f c Fraction of lava that has crystallized f cool Fraction of pyroclasts that reach the surface cooled to ambient f f Fanning friction factor f g Gas volume fraction in magma f hot Fraction of pyroclasts that reach the surface with no cooling f m Lava crystal volume fraction for close packing, 0.45 f o Lava crystal volume fraction for onset of yield strength, 0.021 f r Mass fraction of regolith clasts entrained into lava flow f s Mass fraction of solids in flow after regolith entrainment f trans Total mass fraction of solids in flow at turbulent-to-laminar Specific volume fraction of liquid s Stefan-Boltzmann constant, 5.67 × 10 −8 W m −2 K −4 s h Specific heat of lava, temperature dependent w General value of total lava flow width w i Initial total lava flow width x Distance along flow from vent y Depth below surface z Thickness of thermal boundary layer at top and base of flow ΔE C Energy per unit mass released by crystallization ΔE D Energy per unit mass of lava due to viscous dissipation ΔE G Heat transferred to the ground per unit mass of lava ΔE R Energy loss by radiation per unit mass of lava Δx Increment in travel distance of flow element Λ Heated ground layer depth / flow base thermal boundary layer depth Φ Ratio L h /(T liq − T sol ) α Slope of ground on which lava flows β Bulk density of magma ε Emissivity of lava, 0.96 η Newtonian liquid viscosity η p Plastic viscosity of lava η trans Lava viscosity where turbulent-to-laminar transition occurs κ Thermal diffusivity of lava, 10 −6 m 2 s −1 λ Dimensionless wall film thickness f Gas bubble diameter ρ Magma or lava liquid density σ Surface tension of basalt-water interface, 0.37 J m −2 τ c Time since turbulent-to-laminar transition τ d Exponential flux decay constant τ e Exposure time of lava on turbulent flow surface τ f Time taken by lava to reach a given distance from the vent τ m Time needed to initiate substrate melting at any given distance from the vent τ s Time after eruption start when flow front stops advancing τ str Flight time of pyroclasts in Strombolian phase ORCID iDs Lionel Wilson https:/ /orcid.org/0000-0003-3284-3515