The Control of Lava Rheology on the Formation of Lunar Sinuous Rilles by Substrate Thermal Erosion: Topographic and Morphometric Relationships with Eruption Rates, Erosion Rates, Event Durations, and Erupted Magma Volumes

We describe a model of the fluid mechanics and thermodynamics of turbulent lava flows capable of thermally eroding sinuous rille channels on bodies without atmospheres. The model assumes Bingham plastic rheology for the lava and shows how the effects of radiant cooling and consequent crystallization change the rheology and control the point at which turbulence ceases. A correlation is found between magma volume eruption rate and the length and width of the eroded rille channel. Thus, simple measurements of rille length and width in images can provide reliable estimates of magma eruption rates. The model also predicts rille floor erosion rates, so that if rille depths are measured, eruption durations and erupted magma volumes can also be found. The model is applied in detail to six well-preserved lunar rilles and more generally to a published catalog of 214 lunar rilles. We find that rille-forming eruptions have magma volume eruption rates of a few times 104 m3 s−1, durations of up to 3 months, and erupted magma volumes up to ∼200 km3, consistent with theoretical predictions of basaltic magma ascent and eruption from deep mantle sources. The key requirement for rille formation, rather than mare lava flow deposit formation, is the turbulent eruption of a sufficiently large volume of low-viscosity magma at a sufficiently low eruption rate.


Introduction
We define sinuous rilles (e.g., Figure 1) to be continuously unroofed channels having a monotonic decrease in depth with increasing distance from their sources (Schubert et al. 1970).The spatial correlation of the sources of lunar sinuous rilles with volcanic provinces (Hurwitz et al. 2013) has led to the general acceptance that sinuous rille channels as defined here are the result of dominantly thermal erosion of the preexisting surface by lava flows.Hulme (1973) pointed out that thermal erosion efficiency would be maximized if the lava motion were turbulent.Carr (1974) showed that thermal erosion was still possible, though less efficient, with laminar flow.Hurwitz et al. (2010) discussed the relative importance for typical planetary surface materials of thermal and mechanical erosion by flowing lava, finding that thermal erosion dominated channel formation on slopes much less than 0.1 radians.Slopes for the rilles studied in detail here range from 0.007 to 0.012 radians, so thermal erosion is likely to dominate.In this paper we assume that the lava flows eroding lunar sinuous rille channels left their vents flowing in a turbulent manner.We show that this assumption is consistent with our current understanding of the rates and temperatures at which magmas arrive at the surface of the Moon and with the lengths and depth of sinuous rille channels.
Sinuous rilles are generally easy to distinguish morphologically from partially collapsed lava tubes and lines of vents along the strike of a dike that breached the surface in multiple places (Head & Wilson 2017).Specifically, we require the rille channels in our study to be continuous, and continuously unroofed, for their entire length.For most of their lengths, most sinuous rilles meander.An obvious source of such meandering on the Moon is topographic irregularity caused by ubiquitous impact cratering of the surface.There is evidence that some of the meandering is due to inherent fluid-mechanical instabilities in the motion of the lava producing the channel or spatial variations in the substrate properties (Cataldo et al. 2019).Some sections of rille channels are straight, with evidence of the rille-forming flow having been temporarily captured by a preexisting tectonic depression such as a graben, which was filled with lava and then drained again as lava overflowed and recommenced erosion.Many sinuous rilles have a recognizable source structure-commonly a circular to elliptical depression, more rarely an elongate fissure.Figure 1 shows examples of both source configurations.
The lava in an initially turbulent lava flow erodes its substrate (Figure 2) at a decreasing rate with increasing distance from the vent as the lava cools (Hulme 1973).We find that the rheological changes created by cooling-induced crystallization eventually cause the motion to become laminar, in agreement with the findings of Williams et al. (1998Williams et al. ( , 2000a)).The erosion rate then becomes negligible, and from this point onward a conventional lava flow deposit with positive topography will be formed (Figure 2).If this distal lava reaches a depression and accumulates to form a deepening pond, it is likely that the lava will begin to backfill the eroded channel (Figure 3), so that the depth of the channel, previously decreasing steadily with distance from the source, will abruptly begin to decrease rapidly to zero beyond the point to which it is backfilled.This appears to have happened to both rilles in Figure 1.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Individual lunar mare lava flows with lengths >300 km, widths of ∼40 km, and thicknesses of 10-20 m (Schaber 1973;Robinson et al. 2012) have volumes of ∼200 km 3 .These mafic magmas reached the lunar surface through the anorthositic crust, in which they were commonly negatively buoyant, with minimal physical, thermal, and chemical interaction (Finnila et al. 1994), implying high magma driving pressures, high rise speeds, wide dikes, and near-liquidus eruption temperatures.A mechanism providing these conditions involves the formation of voluminous dikes that grow slowly from the tops of diapiric partial melt zones at depths of several hundred kilometers in the mantle (Wilson & Head 2017, their Figure 1).The dikes disconnect from their sources, migrate upward very rapidly owing to their positive buoyancy in the mantle, and stall centered on the crust-mantle boundary as they reach neutral buoyancy.If such dikes have a great enough vertical extent, their tops penetrate completely through the crust, producing a fissure vent at the surface (Head & Wilson 2017, their Figure 16).Alternatively, a laccolith-like magma body may be formed at the neutral density level at the base of the crust by the accumulation of melt from multiple dikes that had too small a vertical extent to allow them to reach the surface (e.g., Wilson & Head 2023a).Arrival of an unusually voluminous dike could overpressurize the accumulated magma body and drive a dike from its roof to the surface.A similar laccolithic 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, their Figure 2).Based on the buoyancy of mafic magma in its mantle host rocks, Wilson & Head (2017) predicted that dikes growing from mantle diapirs could have vertical extents up to 60-80 km and contain at least 500 km 3 of magma.When the tip of an 80 km tall dike reached the surface while the lower part of the dike was still rising through the mantle, the initial magma volume eruption rate would have been large, ∼10 6 m 3 s −1 .The rate would decrease to ∼10 5 m 3 s −1 over the course of a few days as the dike approaches equilibrium with the density structure of the lithosphere and then decrease further, to ∼10 4 m 3 s −1 , over the course of about a day as buoyancy equilibrium was reached (Wilson & Head 2018).Subsequently, lithospheric forces, compressive during the latter half of lunar geologic history (Solomon & Head 1979), would cause the dike to become narrower, with magma discharge continuing at ∼10 4 m 3 s −1 until no more uncooled magma was available to be erupted.At the other extreme, a dike just able to reach the surface on the lunar nearside where the crust is ∼30 km thick (Wieczorek et al. 2013) would have had a vertical extent a little more than 60 km and a volume of ∼100 km 3 , most of which would eventually be erupted.However, the high driving pressure present at the start of the eruption in the cases of the largest dikes (Wilson & Head 2018) would be absent for these smaller dikes, and most of the eruption would take place at a volume flux of a few times 10 4 m 3 s −1 in an eruption lasting at least several weeks.This prolonged duration is much greater than the time needed for the ground beneath the flow to be heated to its solidus, a few days (Hulme 1982), and guarantees that significant thermal erosion will occur, leading to sinuous rille formation.The above scenario for lunar volcanic eruptions, with large lava volumes being extracted from deep sources and residual magma in the conduit being solidified by the end of the eruption, suggests that multiple eruptions from the same vent were unlikely and therefore rare (Head & Wilson 2017).
In his model of sinuous rille formation by thermal erosion of its substrate by a turbulent lava flow, Hulme (1973) allowed for the temperature variation of the viscosity of the lava but assumed a Newtonian rheology at all temperatures.Hulme (1982) explored models of the cross-sectional shapes of lava flows having Bingham rheology, using these to discuss controls on the widths of the central channels of the lava flows that erode sinuous rilles and hence of the rilles themselves.However, the non-Newtonian rheology was not applied to the flow dynamics during the erosion process, and the implications for the relationship between rille length and eruption conditions were not explored.The lunar rille formation model of Williams et al. (2000a) incorporates the influence of both decreasing temperature and the consequent increase in crystal formation on lava viscosity but does not explore the onset and growth of non-Newtonian behavior.The same is true of analyses of thermal erosion during komatiite eruptions on Earth (Huppert & Sparks 1985;Williams et al. 1998Williams et al. , 2001b) ) and basaltic eruptions on Io (Williams et al. 2000b(Williams et al. , 2001a) ) and Mars (Williams et al. 2005).In modeling the rapid eruption of mare lava flows that are emplaced too quickly for significant thermal erosion to occur, Wilson & Head (2023b) pointed out that the cooling efficiency of a turbulent lava flow is so great that the consequences of increasing crystal content will dominate the rheology of the lava and eventually force the cessation of turbulence.This is a critical issue in the case of rille-forming flows since turbulence is required for significant erosion.We therefore develop a model of the thermal erosion of rille channels by turbulent lavas whose rheology becomes progressively non-Newtonian with distance from the vent.We modify the model of Wilson & Head (2023b) for lower eruption rates, smaller volumes, and longer-duration eruptions, taking detailed account of the thermal interaction between the flow and its substrate.

Fluid Dynamics of Rille-forming Lava Flows
Lunar lavas erupted at near-liquidus temperatures will have Newtonian rheology, i.e., the same viscosity at all shear rates and zero shear rate at zero shear stress.However, the rapid bulk cooling experienced by a turbulent flow produces a significant crystal content, and even small amounts of crystals can cause non-Newtonian behavior (Dobran 2001).We adopt the simplest non-Newtonian rheology, that of a Bingham plastic, possessing a finite yield strength, Y, that must be exceeded before any shear takes place and a constant plastic viscosity, η p , defined as the rate of change of shear stress with strain rate (Skelland 1967).Hulme (1973) found that the viscosity η l of the purely liquid component of lunar lavas could be well approximated as a function of the temperature T by functions of the form As our baseline lunar lava composition, we use the Apollo 12 sample 12002, a low-TiO 2 basalt, which has the advantage that it crystallizes mainly a single component, olivine, over the cooling range relevant to flow modeling (Williams et al. 2000a).This means that the crystal content f c increases uniformly from 0 to 1 as the lava temperature decreases from the liquidus temperature T liq to the solidus temperature T sol .The data given by Williams et al. (2000a) are well approximated by Equation (1) with B = 1477.15K and q = 10.0207,implying a liquidus viscosity of 0.227 Pa s.The liquid viscosity must be multiplied by a correction factor f v , which is a function of the crystal content f c , to give the plastic viscosity η p .We follow Williams et al. (2000a) in using the Roscoe-Einstein equation, for small values of f c (Jeffrey & Acrivos 1976) and the function for larger crystal fractions (Pinkerton & Stevenson 1992).Numerous expressions exist in the literature for the variation of lava yield strength with crystal content (e.g., Gay et al. 1969;Ryerson et al. 1988;Dragoni 1989;Zhou et al. 1995;Saar et al. 2001;Ishibashi & Sato 2010;Mueller et al. 2010;Pantet et al. 2010).We use the following function to approximate the geometric mean of these data: We assume that the lava eroding sinuous rilles leaves the vent flowing in a turbulent fashion.The mean flow speed U is therefore given by where g is the acceleration due to gravity, D is the flow depth, α is the slope of the surface on which the flow moves, and f F is the Fanning friction factor.Data compiled from experimental studies of the flow of Bingham plastics by Skelland (1967) show that f F can be approximated as a function of the Reynolds number, Re, by where the Reynolds number is defined as and ρ is the bulk density of the flowing lava.Turbulence in a flow is maintained as long as the Reynolds number Re is greater than some critical value, Re crit .For Newtonian fluids, Re crit is a constant with a value close to 2200.For Bingham fluids, Re crit is not a constant but instead varies with the viscosity and yield strength.The relationship is defined in terms of the dimensionless Hedström number and by fitting data in Skelland (1967, his

( ) =
This function is valid for He 15,000; at smaller He the rheology rapidly approaches Newtonian and Re crit becomes constant at ∼2200.Note that He depends on the ratio (Y / η p 2 ).Both these parameters increase as the lava cools and its crystal content increases, but in different ways.Predicting the behavior of a cooling and crystallizing turbulent flow is therefore not trivial.In practice, we find that as crystallization begins, both He and Re crit increase rapidly from low values, whereas Re decreases from an initially large value, ∼50 times greater than Re crit .Eventually both He and Re crit decrease slowly, whereas Re decreases rapidly, eventually falling below Re crit and ending turbulence.Equation ( 8) is based on Figure 3.2 in Skelland (1967), which was derived by fitting simple functions to experimental data.Examples of such data in Skellandʼs Figures 5.5 and 6.10 show that the theoretical curves underlying Equation (8) always slightly underestimate the friction factor in a transition region between laminar and turbulent motion that occupies a range of Re spanning a factor of ∼3.Comparing this with the ∼50-fold decrease in Re between the vent and the transition, and noting that Re varies logarithmically with distance from the vent, this implies that the theoretical model overestimates rille lengths, but by no more than 3%.

Thermodynamics of Rille Formation
Hulme (1973) assumed that heat loss from a turbulent lava flow with uniform bulk temperature T would be dominated by radiation from the upper surface of the flow, with surface material continuously replaced by material stirred from below by the turbulence.A key element of this assumption is that the lava does not form a stable cooling crust on its surface.All other treatments of sinuous rille formation (Huppert & Sparks 1985;Williams et al. 1998Williams et al. , 2000aWilliams et al. , 2000bWilliams et al. , 2001bWilliams et al. , 2001aWilliams et al. , 2005) ) have explicitly or tacitly assumed that a stable crust is able to form on turbulent flows, but we see no justification for this when the flow is in a vacuum with no medium (such as an atmosphere or water) in contact with the flow surface.Therefore, like Hulme (1973), we infer that the radiated energy loss per unit mass, ΔE R , as the flow advances a distance Δx is where T a is the average lunar ambient temperature, ∼180 K (Williams et al. 2017), s is the Stefan−Boltzmann constant, 5.67 × 10 −8 W m −2 K −4 , and ε, U, and D are the emissivity, mean speed, and thickness of the lava, respectively.We assume ε = 0.96 (Davies et al. 2011) for the high temperatures involved.
Heat is also lost from the lava owing to the assimilation of the melted substrate.The heat per unit mass, ΔE G , removed in this way is where T mob is the temperature at which the heated rock can be mobilized; (dy/dt) is the rate of advance of the melted ground, i.e., the erosion rate; and c h is the specific heat of the lava in the melting temperature range.Previous treatments generally assumed that T mob was the solidus temperature of the heated substrate.However, rocks melt progressively between the solidus and liquidus, and some finite degree of melting is needed before the partially melted mass can be sheared (Huppert & Sparks 1985).We assume that the substrate rock has similar thermal properties to those of the eroding flow and define T mob to be the temperature at which the corresponding crystal content causes the yield strength of the substrate rock given by Equation (3) to be equal to the shear stress at the base of the eroding flow, (ρgDsinα).The melted substrate is heated to the mobilization temperature T mob , requiring the sensible heat to raise it to that temperature from the ambient temperature plus the appropriate fraction of its latent heat between the solidus and mobilization temperature, so that where h C is the heat transfer coefficient and c c is the average specific heat of the lava between the ambient and solidus temperatures.
Care is needed in specifying h C .All treatments of turbulent heat transfer take the form where K, a and b are dimensionless constants and the Prandtl number Pr is defined as The Nusselt number, Nu, for a fluid flowing in a pipe of diameter D p is defined as (h C D p )/k, where h C is the heat transfer coefficient and k is the lava thermal conductivity.For other geometries the pipe diameter must be replaced by the hydraulic diameter, D h , of the system, defined to be 4 times the cross-sectional area of the moving fluid divided by the length of its wetted perimeter in contact with stationary surroundings.For the lava in the central channel of a flow with depth D and width W between its levees, the cross-sectional area is WD and the wetted perimeter, assuming no rigid crust on the surface of the flow, is W + 2D.D h is therefore equal to (4WD)/(W + 2D).In all cases of interest W is very much greater than D, and so D h is closely equal to 4D, and the Nusselt number becomes  (1973) formula assumed that all components of the system were at similar temperatures and, since there is a strong variation in temperature in the boundary layer where lava is in contact with the ground, chose a formula from Kakac et al. (1987).Unfortunately, they appear to have confused the lava flow depth and hydraulic diameter, so that their formula should be h kD 0.00675 Pr Re , 16 where η b is the viscosity of lava in the boundary layer.Choosing a value for η b is not trivial.At a temperature midway between the lava core temperature and the solidus, η b could be two to three orders of magnitude greater than η p , causing the factor (η p /η b ) 0.14 to be 0.52−0.38.Many of these issues are discussed by Taler & Taler (2017), who propose new equations, validated against experimental data in the literature, for three ranges of Prandtl numbers, 0.1 < Pr < 1, 1 < Pr < 3, and 3 < Pr < 1000, in each case valid for Reynolds numbers in the range 3 × 10 3 to 10 6 .This range of Reynolds numbers overlaps the range relevant to sinuous rille formation, but the rille lavas have Prandtl numbers higher than those studied by Taler & Taler (2017), generally extending between 10 3 and 10 5 .We have used the trends of the coefficients K, a, and b found by Taler & Taler (2017) Figure 4 compares the formulae given by Equations (15), ( 16), and (17) for a range of typical conditions along a turbulent flow as Pr increases and Re decreases with distance from the vent.The lava has thermal conductivity 0.5 W m −1 K −1 , and the flow is 3.7 m deep with a speed of 6.2 m s −1 .The formula adopted here yields higher erosion rates than those of Hulme (1973) and Williams et al. (2000a) near the vent and lower rates near the end of the rille.
In steady flow the kinetic energy of the lava is constant, and so the energy per unit mass produced by viscous dissipation as work is done against friction with the ground, ΔE D , is mainly derived from the potential energy lost by the lava: allows ΔT to be evaluated at each spatial integration step Δx.

Model Implementation
The above equations have been implemented as a spreadsheet program and can be applied to any given rille as follows.A lava composition is inferred based on any available data (spectroscopic; sample) relevant to the rille location.This provides the density, specific and latent heats, and thermal conductivity, assuming that these have been measured.The lava composition can also be used to calculate the crystallization history, which leads to the expressions for the plastic viscosity and yield strength of the lava as a function of its crystallinity and temperature.From Lunar Reconnaissance Orbiter (LRO; Chin et al. 2007) image data (LROC; Robinson et al. 2010) the length, L, of the rille and its width, W, near its source vent can be measured, and LRO Lunar Orbiter Laser Altimeter (LOLA; Smith et al. 2010) topographic data give the slope of the ground on which the rille was formed and the depth of the rille in the immediate vicinity of its source.Some judgment is needed in choosing this depth since the lava pond or fissure that is the source of the flow may show evidence of lava drainback or other modification at the end of the eruption.Additionally, both the depth and the width of the rille will have been modified from their original values by subsequent regolith development and by larger-scale impact cratering.To take account of regolith formation, we measure the rille width at both the top and the base of the rille walls and take the average as the original width of the rille.Where some of the rilleʼs original flat floor still exists, this average supplies a reasonable estimate of the original rille width, and the measured depth to the flat floor provides a good estimate of the original depth.Where the regolith debris slopes meet and no flat floor remains, the measured depth is a lower limit on the original depth and half of the rim-to-rim width is an overestimate of the original width.Finally, regolith development inevitably makes the rille channel difficult to locate once the thickness of regolith that has accumulated since the eruption becomes significantly greater than the channel depth.This acts to reduce the apparent lengths of the oldest rilles.
The inputs to the spreadsheet program, in addition to the material properties just listed, are the eruption temperature and a trial value of the volume eruption rate, F. The eruption temperature is assumed to be close to or a little less than the liquidus, based on the evidence for minimal interaction between rising mafic magmas and the anorthositic lunar crust (Finnila et al. 1994).
The continuity equation for the flow is defined by and dividing the assumed value of F by the measured value of W gives the initial value of the product (UD) at the vent.Equation (6) then gives the Reynolds number for the flow as it left the vent, and Equation (5) gives the corresponding friction factor.Combining Equation (21) with Equation (4), we find and substituting this into Equation (21), so that the mean flow speed and the flow depth can now be found separately.At this point the values of all of the variables characterizing the lava in the flow have been determined, and Equation (20) can be used to find the decrease in the lava temperature as the lava advances by a suitably small integration distance step Δx.This changes the rheological parameters and hence the flow speed and depth, together with the Hedström, Reynolds, and critical Reynolds numbers.For flows tens of kilometers long, Δx = 100 m is adequate.The process is repeated until the Reynolds number has decreased to become equal to the current critical value.

Results
We first give an example for the lower rille in Figure 1, numbered 40 in the catalog of Hurwitz et al. (2013) and 3 in the catalog of Oberbeck et al. (1971).We have measured its length, near-vent width and depth, and average substrate slope using LRO LROC images and LOLA topography accessed via Quickmap,3 the values being 57.4 km, 807 m, 209 m, and 0.012 1 radians.We implement the model using the properties of the low-TiO 2 mare basalt determined by Williams et al. (2000a).The dense rock equivalent volume eruption rate required to cause erosion of the rille channel to stop at the measured 57.4 km from the vent is F = 1.36 × 10 4 m 3 s −1 .The initial mean flow speed and flow depth are 5.8 m s −1 and 2.9 m, respectively, and the initial Reynolds number is 4.60 × 10 5 .The eight panels of Figure 5 show the change with distance from the vent of the parameters characterizing the eroding flow.Initially the yield strength is negligible, the Hedström number is zero, the critical Reynolds number is 2200, and the lava is fully turbulent and Newtonian.As the lava cools and its crystal content increases, non-Newtonian behavior begins 1.3 km from the vent with a rapid increase in He and Re crit along the flow.After rising to maxima at ∼20 km from the vent, both He and Re crit start to decrease slowly.In contrast, the Reynolds number Re has been decreasing steadily along the flow, and its value converges with that of Re crit at 57.4 km.At this point, turbulence ceases.The crystal content is then 45.2%, and the yield strength and plastic viscosity of the lava are 105 Pa and 24.1 Pa s, respectively.The values near the vent of the terms ΔE D /Δx, ΔE R /Δx, and ΔE G /Δx for this flow are 0.019 6, 9.52, and 0.287 J m −1 , respectively.The term ΔE C /ΔT is 3660 J K −1 , so that over an integration step length of Δx = 100 m the temperature decreases by 0.267 K. Thus, ∼97% of the heat loss is by radiation from the flow surface, and ∼3% is by loss to the underlying surface.The erosion rate of the preexisting surface beneath the flow near the vent is 22.7 μm s −1 , i.e., almost exactly 2 m day −1 .Dividing the measured 209 m rille depth by this rate implies an eruption duration of 107 days; multiplying this duration by the 1.36 × 10 4 m 3 s −1 eruption rate implies an erupted volume of 126 km 3 .
The above results change if a different lava composition is assumed.Using the properties of the high-TiO 2 Apollo 17 basalt also reported by Williams et al. (2000a), the eruption rate is reduced by ∼12% to 1.21 × 10 4 m 3 s −1 .The erosion rate increases by 17% to 27 μm s −1 (2.3 m day −1 ), the eruption duration decreases by 16% to ∼91 days, and the erupted volume decreases by 28% to 95 km 3 .
We have used the low-TiO 2 composition data to model in detail five additional rilles.The results are shown in Table 1.With one exception, eruption rates are of order a few × 10 4 m 3 s −1 , erosion rates are all in the range 17-23 μm s −1 (1.5-2.0 m day −1 ), eruption durations are 9-15 weeks, and erupted volumes of lava needed to erode the rile channels are 17-203 km 3 .These results are entirely consistent with the predictions of the Wilson & Head (2017) model of magma supply from the lunar mantle in dikes just vertically extensive enough to penetrate completely through the lunar crust.
It is interesting to compare the volumes of lava needed to erode the rille channels with the volumes of bedrock removed to create the channels.Table 2 approximates the missing rock volume by assuming that the rille width is approximately constant.The erosion rates in panel (h) of Figure 5 imply that the average rille depth is expected to be ∼25% of the depth next to the vent, so the eroded volume is assumed to be onequarter of the product of the rille length with the depth and width near the vent.Table 2 then shows that the amount, C, by which the lava flowing beyond the end of the rille channel has been contaminated by the older rock amounts to 1%-2% by volume.
Figure 6 shows the depth profile of the second rille in Table 1, which has the smallest fitted eruption rate.Also shown are two scaled versions of the erosion rate curve-the equivalent of panel (h) of Figure 5 for this rille.The erratic scatter of the 11 depth measurements illustrates the problems of measuring rille depths and extrapolating the depth to zero to find the true length of the eroded rille.This is particularly difficult when the distal part of the channel has been backfilled by ponded lava.As is common for rilles, in places this rille is unusually shallow where debris from nearby impacts has fallen into it.In other places the rille is anomalously deep where a crater has formed on the rille floor.Additional problems include deformation of the rille by post-formation wrinkle ridge development.Finally, the shallowest, distal parts of the rille will be most influenced by the general lateral regolith mixing caused by impact cratering.The measurements shown in Figure 6 were made at a mixture of locations chosen to either avoid or deliberately include these problems.The most likely immediately post-eruption profile is that represented by the trace of the deepest parts of the rille.We have therefore taken the theoretical erosion profile for this rille and expanded its xaxis by various factors to look for a good fit with the deepest parts of the rille.Figure 6 includes two of these expanded profiles, the expansion factors being 1.3 and 1.7.Comparing the observed locally deepest points with the theoretical curves suggests that a factor of ∼1.6 may be appropriate, implying that the original rille length was ∼1.6 times the measured 32.8 km length of the rille, i.e., ∼52.5 km.The erupted volume flux needed to produce this rille length is 5.20 × 10 3 m 3 s −1 , rather than the original 3.25 × 10 3 m 3 s −1 estimate, also larger by a factor of 1.6.
An important correlation emerges from Table 3, which reproduces the measured lengths and widths of the rilles and the deduced eruption rates from Table 1.Also given are the results of dividing the eruption rates by the products of the lengths and widths of the rilles, now both given in km.The resulting numbers are extremely similar: the average is 301 ± 9 m 3 s −1 km −2 .The implication is that if the width and length of a rille channel are measured, bearing in mind the earlier caveats about estimating the extent of any drowned distal part of the rille, the eruption rate of the lava eroding the rille can be found from where F is measured in m 3 s −1 and W and L are measured in km.This result holds specifically for lava with the composition of the low-TiO 2 Apollo 12 sample 12002.We have produced the equivalent of Table 3 using data for the high-TiO 2 composition Apollo 17 sample 74220.In this case the resulting correlation is Note.Hurwitz number and Oberbeck number are the numbers assigned to the rilles in the catalogs of Hurwitz et al. (2013) and Oberbeck et al. (1971), respectively.Lat. and Long.are the latitude and longitude of the center of the vent, respectively; L is the visible rille length measured along the thalweg; W and D are the width and depth, respectively, of the channel at the exit from the vent structure; α is the slope of the pre-eruption surface measured along the thalweg; F is the volume eruption rate of magma; dydt is the erosion rate of the ground beneath the flow; τ is the implied eruption duration; and V e is the total volume of magma erupted.Note.Hurwitz number: see Table 1.V e is the volume of erupted lava that eroded the rille channel.V r is the volume of pre-eruption surface rock eroded during the eruption.The ratio of these volumes implies the amount C by which the lava flowing beyond the end of the rilles channel has been contaminated by the older bedrock.Again, the composition of the lava cannot be ignored in modeling rille formation.Even so, we suggest that if all that is known about a rille channel is its length and width, then a much better than order-of-magnitude estimate of the lava eruption rate during the event that formed it can be found from using the average of the above two equations, F = 286 WL.

Discussion
The above findings allow us to capitalize on the data collected by Hurwitz et al. (2013), which include measurements of the lengths and average widths of 214 lunar rilles.Figures 7(a)-(c) are histograms of the volume eruption rates, eruption durations, and erupted volumes of lava deduced from rille lengths, widths, and depths for all these rilles except the most extreme examples, discussed below.To allow for the difficulty of locating the termini of the rilles, a factor of 1.4 has been applied to the average of the eruption rates derived from Equations (24) and (25).A few adjustments have been made to the catalog measurements where newer data, especially for Rimae Sharp and Mairan, were available from Qian et al. (2021).All the inferred eruption rates are consistent with the eruption model presented here.Even the smallest value of F, 300 m 3 s −1 , derived from the 3 km long rille number 158 in the Hurwitz et al. (2013) catalog, implies that the lava forming this rille had a Reynolds number of 10,500 as it left its source depression, well above the 2200 threshold for turbulence provided that its temperature was close to its liquidus.However, it is clear from examination of images that many of the short rilles have been partly buried by later eruptive or impact activity or may have roofed-over sections.Recall that our model is only intended to be applicable to channels that are unroofed over their entire length.In addition, 131 rilles do not show any clear evidence of a vent or source depression and so must be partly obscured.These observations imply that the histogram peaks in Figure 7 for eruption rates less than 10 4 m 3 s −1 , durations less than ∼30 days, and volumes less than ∼30 km 3 likely overestimate the proportion of eruptions that had such low values.
Rima Schroeter is not included in Figure 7.Its length of 176 km and near-vent width of 4.3 km imply an eruption rate of 3.05 × 10 5 m 3 s −1 .These properties, as well as those of the four rilles with the next-largest inferred eruption rates, are also listed in Table 4.As is the case for the rilles requiring smaller lava eruption rates, parameters for all these rilles are consistent with the model developed here, although these examples are predicted to have eruption rates more comparable to those of mare lava flows than sinuous rilles.This implies that there should be very extensive "regular" lava flow deposits connected to the distal ends of these rilles.It is beyond the scope of this analysis to assess evidence for these, but careful mapping might reveal their presence.An example of one such deposit was identified for Rima Prinz by Hurwitz et al. (2012).
The most striking features of Table 4 are the very large volumes predicted for the flows eroding Rima Schroeter and Rima Plato.The value for Rima Plato is at the extreme upper limit of the range of single-eruption magma volumes predicted by the Wilson & Head (2017) mantle dike model, and the value for Rima Schroeter is ∼4.5 times larger than that limit.These extreme volume predictions arise because of the unusually great depths of the rille channels.Rima Plato is located on the ejecta blanket of the Plato impact crater, and it is possible that mechanical erosion of the poorly consolidated ejecta deposit may have played a part in causing the great depth of this rille, as is commonly seen in short, deep sinuous rilles in the Orientale basin that have sources in nonmare material (Whitten et al. 2011).Rima Schroeter is unusual in many respects (e.g., Jawin et al. 2016) and clearly requires a more detailed study of its unusual morphology (e.g., Wilson et al. 2021).

Conclusions
1. We present a model for the thermal erosion of sinuous rille channels by turbulently flowing lava.The lava loses heat by radiation from its constantly stirred surface, and no continuous insulating crust is formed.The model emphasizes the importance of the evolving non-  Newtonian rheology of the lava as it loses heat and crystallizes.2. The changing rheology provides the specific reason for the cessation of turbulence and hence the end of significant substrate erosion and allows the length of the rille channel to be linked to the volume eruption rate of lava and the width of the channel.Readily available rille length and width measurements from images therefore allow estimates of the lava eruption rates required to form individual rilles.High volume eruption rates of lava are implied, and these compensate for the efficient heat loss from the turbulent lava surface and allow rilles of the observed lengths to be formed without the need for insulating crusts on the flows.3. The model also predicts substrate erosion rates, so that topographic measurements of rille channel depths from laser altimetry, stereometry, or shadow length measurements can provide estimates of eruption durations.
Multiplying the duration by the volume eruption rate then provides the volume of lava erupted and emplaced as a conventional lava flow deposit beyond the end of the rille channel.4. Volume eruption rates, eruption durations, and erupted lava volumes have been derived for the flows eroding the 214 rilles whose morphological properties were measured by Hurwitz et al. (2013).For rilles that have survived major modification by later eruptions, typical eruption rates are a few times 10 4 m 3 s −1 , durations range from a few weeks to 6 months, and volumes range from several tens to several hundreds of km 3 .These values are consistent with a model (Wilson & Head 2017) for the ascent and eruption of lunar magmas derived from deep mantle sources. 5.The key difference between an eruption that emplaces one of the "classic" mare basin lava flow deposits like those in Mare Imbrium (Schubert et al. 1970) and an eruption that leads to the formation of a sinuous rille is the arrival at the surface of a mantle-derived dike erupting a sufficiently large volume of magma at a sufficiently low eruption rate that the eruption continues for at least a few weeks.The mare lava flow eruptions have ended before the substrate has been heated to the point of efficient melting, whereas the rille-forming eruptions continue long enough that most of the eruption takes place in the efficient substrate melting regime.

Figure 2 .
Figure 2. Turbulent flow leaves vent and transitions to laminar after radiative cooling causes increase in crystal content and increasingly less Newtonian rheology.Substrate erosion rate decreases with distance from vent and ceases at turbulent-laminar transition.

Figure 3 .
Figure 3. Cross section and plan view of turbulent flow eroding substrate to form a sinuous rille channel.Lava ponds in depression and may backfill the rille channel, causing onset of rapid decrease in channel depth and leading to underestimation of channel length in images.
the values A = 0.044 56 Pa, f o = 0.022 09, f m = 0.6, and p = 1.773 75.Yield strengths of a few hundred pascals at moderate solids contents are implied, consistent with measurements made on terrestrial mafic lavas by Robert et al. (2014) and Sehlke et al. (2014).

Figure
The parameters K, a, and b in Equation (12) are obtained from experimental literature.Hulme (1973) adopted K = 0.023, a = 0.4, and b = 0.8Williams et al. (2000a) noted that Hulmeʼs sensible and latent heat energy per unit mass released as the temperature decreases by ΔT during c h and L h are the specific heat and latent heat of the lava, respectively.Since ΔE D is a heat gain, whereas ΔE R and ΔE G are heat losses, we have ΔE C = ΔE D − ΔE R − ΔE G , and so balancing gains and losses,

Figure 4 .
Figure 4. (a) Variation of Prandtl and Reynolds numbers with distance from vent.(b) Corresponding variations of the heat flux transferred to the ground using the heat transfer coefficients of Hulme (1973), Williams et al. (2000a), and this work.

Figure 5 .
Figure5.Variation with distance from the vent of eight parameters characterizing the lava flow that eroded the channel of rille number 40 in Table1.

Figure 6 .
Figure 6.Eleven depth measurements along rille number 40 listed in Table1compared with theoretical curves assuming that the measured length has been underestimated by factors of 1.3 and 1.7.
to extrapolate them to the higher Pr values, finding

Table 1
Thermal Erosion Model Results for Six Rilles Note.See Table1for parameter definitions.

Table 4
Properties of Rilles Formed by the Highest Eruption Rate Flows Note. see Table1for parameter definitions.
Stefan-Boltzmann constant, 5.67 × 10 −8 W m −2 K −4 t time x distance along flow from vent y depth below surface ΔE C energy per unit mass released by crystallization ΔE D energy per unit mass pf 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 ΔT change in lava temperature Δx increment in travel distance of flow element α slope of ground on which lava flows ε emissivity of lava, 0.96 η l