## Abstract

Carbonaceous asteroids may have been the precursors to the terrestrial planets, yet despite their importance, numerous attempts to model their early solar system geological history have not converged on a solution. The assumption has been that hydrothermal alteration was occurring in rocky asteroids with material properties similar to meteorites. However, these bodies would have accreted as a high-porosity aggregate of igneous clasts (chondrules) and fine-grained primordial dust, with ice filling much of the pore space. Short-lived radionuclides melted the ice, and aqueous alteration of anhydrous minerals followed. However, at the moment when the ice melted, no geological process had acted to lithify this material. It would have been a mud, rather than a rock. We tested the effect of removing the assumption of lithification. We find that if the body accretes unsorted chondrules, then large-scale mud convection is capable of producing a size-sorted chondrule population (if the body accretes an aerodynamically sorted chondrule population, then no further sorting occurs). Mud convection both moderates internal temperature and reduces variation in temperature throughout the object. As the system is thoroughly mixed, soluble elements are not fractionated, preserving primitive chemistry. Isotopic and redox heterogeneity in secondary phases over short length scales is expected, as individual particles experience a range of temperature and water-rock histories until they are brought together in their final configuration at the end of convection. These results are consistent with observations from aqueously altered meteorites (CI and CM chondrites) and spectra of primitive asteroids. The “mudball” model appears to be a general solution: Bodies spanning a ×1000 mass range show similar behavior.

## INTRODUCTION

Primordial asteroids accreted as a mixture of nebular fines, ice, and coarse spherules (chondrules). Carbonaceous chondrite meteorites sample these asteroids. A subset of these meteorites shows evidence of significant hydrothermal alteration, apparently occurring within their asteroid parent bodies [for example, see the study of Brearley (*1*)]. Counterintuitively, the most chemically pristine meteorites (CI and CM chondrites) have experienced the most pervasive aqueous alteration. CI and CM meteorites contain abundant water and organics and have chemical compositions that are a close match to the solar photosphere—in the case of CI chondrites, at least 40 elements within 10% of the photosphere value (*2*, *3*). The solar/meteoritic ratio is independent of volatility or geochemical affinity. When considering the hydrothermal history of a geological body, it is also independent of solubility. Highly soluble elements, such as Na, K, and F, are also within 10% of solar (*2*, *3*). The unusual chemistry of CI and CM chondrites, abundant water, and mix of complex organics have led many to identify these meteorites as plausible terrestrial planet precursors or candidates for delivery of terrestrial volatiles. A desire to understand their geological history, and the early evolution of water and organics in the solar system, has been a driver for numerous asteroid-scale numerical modeling studies (*4*–*10*).

Hydrothermal alteration transformed the mineralogy of these rocks—CI chondrites are composed of >95% secondary phases, with only a trace of their original anhydrous mineralogy remaining (*11*)—without changing their chemistry. This “isochemical” alteration may result from negligible fluid flow—a closed system at scales larger than 1 mm. However, earlier numerical studies have consistently predicted large-scale (tens of kilometers) water transport in model asteroids (*4*–*10*). This type of convective heat transfer is an effective way of moderating temperatures within an asteroid parent body [for example, see the study of Travis and Schubert (*8*)], something that is particularly important in the case of CI and CM chondrites, because almost all of these rocks appear to have experienced alteration over a narrow temperature range and at relatively low temperatures (<150°C). However, in addition to asteroid-scale convection, which might fractionate aqueous species, a problem with numerical studies is that they often produce only narrow bands within model asteroids that match CM chondrite conditions (*4*–*6*, *8*, *10*). A consideration of the material properties of chondrite precursors (note, not extant meteorites, because hydrothermal alteration and impact processing have—in almost all cases—profoundly changed their material properties) suggests very low permeability (*12*): Earlier geophysical models may have used permeability estimates around six orders of magnitude higher (*12*). As a result, liquid water flow would be reduced to millimeter-scale distances at most, even in a high-porosity, water-saturated asteroid. Low permeability provides a mechanism for isochemical alteration, but in the absence of convective motion of water (a “closed-system” scenario), we are required to posit small (less than tens of kilometers in diameter) primordial parent bodies to moderate temperatures during alteration (*12*). In addition, carbonate grains in some CMs exhibit a large range in δ^{18}O [for example, see the study of Jenniskens *et al*. (*13*)], requiring varying temperatures and water/rock (WR) ratio over short length scales—difficult to achieve in a low-permeability, closed-system environment. This is interpreted by some as suggesting open-system WR interaction and an asteroid that is zoned with respect to O-isotope composition (*5*). Achieving this O-isotope heterogeneity at the thin section scale would require a dense network of channels [an additional benefit would be that fractures and channels would increase permeability and permit efficient heat transfer (*14*)]. However, this type of fracture/channel network feature has not been observed in any chondrite, on any scale). Reaction modeling (*15*) predicts abundant CH_{4(g)} and H_{2(g)} as a consequence of asteroidal aqueous processing, requiring an open system at least with respect to gaseous species. Finally, young CM chondrite dolomite ages [3.4 to 4.5 My (million years) (*16*) after calcium- and aluminum-rich inclusion (CAI) (*17*)] raise questions as to whether a late-accreting body, containing substantial water ice, could even melt and achieve temperatures consistent with dolomite formation (defining *T*_{peak} for most CMs at ~150°C).

It is apparent that a variety of well-supported model constraints, data and observations, and existing numerical approaches are not converging on a holistic scenario for the hydrothermal evolution of primitive, water-rich (CM and CI chondrite) asteroids. This is unfortunate, because a well-founded geophysical model for the aqueous and thermal alteration that occurred in these objects would give us a framework for understanding the evolution of water and organics in the early solar system and the geological history of objects that may have been the precursors to the terrestrial planets. Here, we propose an alternative that reconciles these competing lines of evidence. All previous studies of water/rock interaction in asteroids have assumed that the object is lithified: An anhydrous coherent rock reacts with water. However, there is no a priori reason why this should be the case. On accretion, nebular fines, ice, and chondrules would not be lithified. On melting, chondrules might become a suspended load in a viscous mud. Water and solids would move together. The notion of open or closed system would no longer apply.

## RESULTS

Here, we explore this concept with the MAGHNUM numerical simulator, previously used to model thermal evolution of consolidated carbonaceous chondrite parent bodies (*8*, *18*). New features include particle transport using Stokes-Rubey settling (*19*, *20*), minimum porosity from packing models, mud viscosity, serpentinization, and dehydration. The code tracks chondrules. Dimensionless tracer particles scattered throughout the mud allow us to track movement, temperature, and WR ratio path of dimensionless mud “grains” (which might translate to individual phyllosilicate or carbonate grains in CM chondrite matrix), as they evolve within the body over time. Another novel code feature concerns porosity and permeability. In an unconsolidated particulate bed, especially at low *G*, these parameters are dynamic. Nonuniform flow can generate feedbacks, increasing (or decreasing) porosity and permeability through particulate movement, which in turn can increase (or decrease) flow, eventually creating structure in initially uniform media. The code handles this nonlinear dynamic interplay of porosity and permeability and flow in an unconsolidated bed, at each time step at each node. Model asteroids are unconsolidated mixtures of coarse particles (chondrules) and mud (initially a uniform mix of fines and ice), with starting compositions chosen to provide a comparison to CM chondrite parent bodies at plausible accretion times. We vary mud/chondrule ratio and mud viscosity. Chondrule packing density is the mud-filled porosity between clasts. We typically use a packing density value of 50% (that is, 50% chondrules, with interstices having a volume of 50%, filled by mud), which approximates a terrestrial sand, although we do vary this in some simulations. The viscosity of the mud in our model follows a power law (*21*). Shoji and Kurita (*22*) used an almost identical relationship. Mud volcano viscosities (*23*) lie in the same range, and data from the geologic literature (*24*, *25*) bracket the other (*21*, *22*) curves. The viscosity of mud shows considerable dependence on particle type and fines content/load fraction [for example, see the studies of Zoporowski and Miller (*23*) (and references therein), Fritz and Moore (*24*), and Grim (*25*)]. To illustrate the impact of a higher mud viscosity, we ran cases with a 10-fold increase in viscosity. We still see convective motion, reduced somewhat in strength but qualitatively the same. We also vary WR ratio [a bulk initial WR of 0.6 and 1.0, bracketing estimates for CM chondrites (*26*)]. WR varies over time and affects mud viscosity, as water is consumed in alteration reactions. At 60:40 matrix/chondrule mix, a bulk initial WR of 0.6 and 1.0 translates to an initial matrix WR of 1 and 1.67. Simulations span asteroid radii from 50 to 500 km. Although the canonical view is that chondrite parent bodies accreted a population of chondrules that had already been sorted in the disk (and sorted chondrule populations in chondrite groups that have low matrix abundance testify that this occurred), to explore whether mud convection is also capable of segregating chondrules based on size, we begin with an unsorted population made up of 1-mm, 0.5-mm, 0.1-mm, 50-μm, and 10-μm particles in a 10:20:30:30:10 ratio. Simulations beginning with a sorted chondrule population produce qualitatively similar results (in terms of *T*_{peak} and thermal structure within the body).

In simulations of unconsolidated asteroids that have a high mud/chondrule ratio, chondrules settle rapidly to form an unconsolidated chondrule + mud “core” and mud ocean (85 to 98% mud) mantle capped with an undifferentiated frozen crust (see the Supplementary Materials regarding the stability of the crust). Because there are no empirical constraints on the initial, bulk, accreted ratio of chondrules to matrix in primordial asteroids (the tacit assumption is that the object was lithified: the observed ratio in the meteorite is taken to be representative of the ratio in the parent body), this is a free parameter. In most simulations presented here, we chose a ratio that allows us to visualize and explore an outcome that is unique to the mud model: differentiation between a core that is richer in chondrules and a mud mantle. However, note that this is not predictive of CM parent body structure. Different choices for initial chondrule fraction and packing density translate to differences in core/mantle volume and core chondrule/matrix ratio. A larger accreted chondrule fraction, lower packing density, and less compaction with depth lead to a reduced—or nonexistent—mud mantle, as illustrated in Fig. 1 in the high viscosity example.

Mud convection can provide an efficient chondrule sorting mechanism (Fig. 1), although, once again, if the object accretes a sorted population via some aerodynamic mechanism [for example, turbulent concentration in the disk (*27*)]—the canonical view—then additional sorting in the parent body cannot occur. There are a number of model inputs that affect the degree of sorting: asteroid radius (smaller objects, with lower gravity, show less efficient sorting), initial viscosity (higher viscosity reduces or removes sorting), WR ratio (lower WR ratio translates to higher viscosity, and because WR evolves over time, so does viscosity), packing density (lower packing density reduces sorting), and accretion time (objects that accrete earlier show stronger convection and more efficient sorting). If sorting occurs (in high-viscosity scenarios, it is absent), then we typically see an inner core region dominated by 1- and 0.5-mm chondrules, surrounded by 0.1-mm particles, with 50 μm concentrated at the edge of the core. Particles (10 μm) are distributed evenly either throughout the object or at the core mantle boundary (Fig. 1). In a number of scenarios, we see particle oscillations, which reflect a feedback between fluid flow and particle transport. Sufficiently strong convective plumes loft coarse particles, which spread out in space (in a manner similar to a chromatograph) because of their size (Stokes-like flow, which depends on particle size). The smaller particles are more susceptible to lofting. Upwelling convective flowing mud lofts particles upward and sideways, with convergence in downwellings. This motion tends to cause particles to be depleted in upwelling regions and concentrated relative to the average in downwelling regions. This type of clustering occurs primarily while particles are falling from the melt front rather than after core formation. It is at this point in the hydrothermal evolution of the parent body when convection is most intense, driven by serpentinization’s exothermic reaction and decay of short-lived radionuclides. Core formation is not instantaneous, taking from 25,000 to >100,000 years, depending on the fine/coarse (f/c) ratio, object size, and WR value. Once coarse particles reach the maximum packing density in any core bed node, the bed becomes immobile (but is not consolidated) at that location.

Convection is not restricted to the mud mantle but extends into the core, where convection consists of mud moving through the porous bed formed by coarse particles. If the accreted f/c ratio (and packing density) allows for differentiation, then we may see plumes in the mud ocean driven partly by plumes in the core that exit to the ocean bottom, as well as cold downwellings dropping off the frozen crust. We can visualize convection by tracking the history of individual dimensionless tracer particles scattered throughout the mud. Considering a scenario that has a 100-km-radius object accreting at 3 My after CAI with a bulk WR of 1.0, Fig. 2 shows individual trajectories for tracer particles that begin at a level of 20 km (Fig. 2A), 50 km (Fig. 2B), and 80 km (Fig. 2C) in this object. Ice melting and serpentinization begin at around 0.9 My after accretion. The figure shows how convection strengthens over time. By 2 My, discrete trajectories for particles that started at all three levels span the entire object—from the core up to the base of the ice lid. Objects accreting after 4 My do not melt. Asteroids accreting at 3 and 3.5 My after CAI take 0.3 to 0.6 My to warm to 0°C from decay of short-lived radioisotopes, where temperature buffers at 0°C for 0.5 to 2 My, whereafter internal heating overcomes the latent heat of ice melting (higher WR and later accretion delay melting). Temperature then increases rapidly, driven by serpentinization. The *T*_{time} trajectory varies depending on the size of the object, but interior to the ice shell, the thermal stratigraphy—the *T*_{avg} profile, where temperature is averaged across all longitudinal coordinates at the same level in the object—is essentially the same, despite the fact that the model asteroids span a ×1000 range in mass (Fig. 3). At *T*_{peak}, a Ceres-class object with WR = 1.0 accreting at CAI + 3 My, a 100-km object with WR = 0.6 accreting at CAI + 3.5 My, and a 100-km object with WR = 1.0 all have an approximately flat *T*_{avg} profile from the interior to the ice shell and similar temperatures of *T*_{avg} = 125 ± 20°C throughout the bulk of the object. Mud convection moderates temperature and minimizes any internal thermal gradient across all model asteroid sizes. Cooling rate is a significant variable with respect to object size: 150 My to reach 50°C in a Ceres-class object (mud convection is still occurring at this time) versus 5 to 10 My in smaller objects. Scenarios with lower f/c ratio, lower packing density, and higher viscosity (Fig. 3E) show somewhat lower temperatures. Mud convection still occurs in these cases, and the *T*_{avg} profile is still flat, but chondrules are not sorted (Fig. 1) and no clear core forms. Peak temperature is cooler than when a core forms, because heat generation is roughly homogeneous throughout the body: Approximately half of the heat generation is in the outer 20% radially; being closer to the surface, the object cools faster.

## DISCUSSION

Our simulations are not intended to be predictive of CM parent body structure (varying f/c ratio, packing density, viscosity, degree of sorting in the accreted chondrule population, or having a sorted chondrule population gives rise to a range of parent body stratigraphies). Some of these variables (for example, accreted f/c ratio) cannot be constrained. However, a robust result is that over a wide range of input parameters, mud convection occurs, generating simulations that are a close match to the conditions for CM chondrite alteration and chronology, from petrographic and isotopic studies. Young CM carbonate ages [formation at ~4 My after CAI (*16*, *28*)] are expected. A planetesimal can accrete earlier, incorporating sufficient short-lived radionuclides to melt, but it will experience an ~1- to 2.4-My delay between accretion and *T*_{max}. Compared to other chondrite groups, CM chondrites record hydrothermal alteration over a relatively narrow (and low) temperature range (*29*). Inferred alteration temperatures from isotopic analyses of CM carbonates are 20° to 71°C (*30*) and 0° to 130°C (*31*). The presence of dolomite indicates a (relatively) high initial parent body temperature, in excess of 120°C (*32*). Reaction modeling of phyllosilicate formation (*15*, *33*) predicts temperatures at the lower end of this range. These observations suggest a CM parent body where internal temperature was moderated and that did not have a steep thermal gradient (not an onion shell). This thermal structure is observed in all of our simulations, and the fact that we see similar *T*_{avg} temperature profiles (and parent body histories) across all simulations suggests that the mudball scenario offers a general model for CM chondrite alteration. CM-like rocks can be extracted from asteroids of a range of sizes that accrete at varying times, with varying WR ratios and viscosities, and from a range of depths within the parent body. This prediction of asteroid-scale homogeneity is consistent with recent observations of CM-like main belt (Ch/Cgh spectral type) asteroids (*34*) and contrasts with highly stratified model objects from previous work (*4*–*6*, *8*, *10*).

Possible asteroid-scale homogeneity gives way to fine-scale chemical and isotopic heterogeneity in CM and CI chondrites. The mudball model provides an explanation for this. Heterogeneity in carbonate oxygen requires formation at varying temperatures and WR ratios, previously interpreted as consistent with large-scale open-system circulation of fluids through the body (*5*). Carbonate heterogeneity is not limited to O-isotopes. “Clumped” isotope and C-isotope data (*30*, *31*, *35*), and chemical data (*35*), also indicate formation over a range of temperatures [even in the same meteorite (*30*)] and in different alteration environments under different physicochemical conditions (for example, redox states)—this heterogeneity is observed at fine scales [hundreds of micrometers (*35*)]. This is surprising if we assume a lithified rock, where neighboring components could be expected to witness similar conditions. However, it is a necessary outcome of a model where grains forming in a range of environments, over extended time scales, are mixed via mud convection. In the mud model grains forming at different times, temperatures and WR ratios will be brought together in their final configuration only at the end of convection. To constrain the degree of heterogeneity, we track the variety of dimensionless tracer particle histories for temperature and WR trajectories for fines, for particles that begin at the same level in an asteroid (Fig. 4). Carbonate formation appears to have been relatively inefficient in CM parent bodies (*31*), occurring throughout the active phase of aqueous alteration. Although this process (continual secondary mineral formation and partial back-reaction) is not explicitly modeled in our simulations, we can visualize its effect using the tracer particle tracking approach (again, we take the example of a 100-km-radius object accreting at CAI + 3 My with bulk WR = 1.0). Average temperatures throughout the object are moderated by convection (Fig. 3). This is observed on a particle basis as well (Fig. 4). By tracking individual tracer particles starting at 60 km (Fig. 4A), we see that in the first 0.5 My after melting, most experience *T*-*t* trajectories in the 120 ± 20°C range. However, a significant minority experience *T*-*t* paths far outside this range—excursions up to 200°C, and 0°C for material eroded from the ice lid. In addition, temperature variability is much greater for particles that begin at higher levels in the object (Fig. 4B, 90 km). By 2 My, this variability has reduced—no particles have trajectories outside ~130 ± 30°C. Serpentinization occurs rapidly (<1000 years), and with it the main modification in WR (Fig. 4, C and D). However, as mud convection strengthens, strong upwelling occasionally degrades portions of the ice lid, leading to transient spikes in WR at high levels in the model asteroid [Fig. 4, C (70 km) and D (90 km)]. In our 100-km-radius example, this period lasts from around 1 to 1.3 My after accretion. It is apparent that there is a significant range in possible *T*-*t* and WR-*t* paths. In addition, continual secondary mineral formation will record changing conditions as the body cools. Compositional and isotopic heterogeneity over short length scales at the end of convection is expected.

There is petrographic evidence in favor of the mudball model. The edges of carbonate grains are abraded in CM chondrites, suggested to be a consequence of fluidization of the matrix (*28*). Chondrule abrasion and layers of different particle sizes and/or particle types have also been described (*36*). Textures interpreted as arising from compaction following impact (*37*) are entirely consistent with mud flow, and a ubiquitous feature of CMs is also explained by it. Rims of fine-grained matrix material surround chondrules in CMs (and other carbonaceous chondrites). Often referred to as “accretionary rims,” they may have acted to cushion impacts, enabling the aggregation of the first solids (*38*). Evidence in favor of a nebula accretionary origin is the relationship between the radius of the clast and the thickness of the fine-grained rim (*39*). The mud model may offer an alternative. CM chondrules will be continually moving through mud during settling or while lofted by strong plumes. Mud convection extends into the chondrule-rich core—mud will be flowing through the porous chondrule-rich bed until convection ceases. In this environment, rim textures may be a meteoritic equivalent of pisoliths in terrestrial sedimentary settings. Chondrule rims in CMs would still be accretionary, but they would have accreted in the parent body rather than the nebula. In addition to explaining this feature of CM petrography, it would have the effect of increasing the matrix/chondrule ratio from the base level, defined by our choice of matrix/chondrule packing density (50%). Whether or not it is a complete rim, any accretion of matrix material onto chondrules will have this effect. It is noteworthy that CM2 chondrites have highly variable matrix abundances, from 58% (Murray) to 85% (Nogoya).

Finally, the mudball model circumvents the isochemical alteration paradox, where the most chemically pristine meteorites are the most aqueously altered. Solubility-related chemical fractionation is only an issue where water is moving through a lithified rock matrix. In a viscous mud, fines and water move together. In all scenarios, convection of mud—whether in a mud mantle or through a porous bed formed by coarse particles—is occurring throughout the entire object up to the ice lid. The system is thoroughly mixed, so alteration will necessarily be isochemical: Local solubility-related fractionation will be reduced or removed by mud flow.

## MATERIALS AND METHODS

### Description of MAGHNUM numerical model

The computational model MAGHNUM combines modeling features described in previous work (*8*, *10*, *18*, *40*), with additional features for particulates flowing in mud and forming an unconsolidated porous, permeable bed. It solves equations for mass conservation, momentum conservation, energy conservation, and solute conservation, with equation of state and constitutive relations, as summarized below. In addition, several issues of interest are discussed.

Mass conservation of H_{2}O. Rate of change of mass of H_{2}O in a volume depends on advection of water(1)where *s*_{ice} is the fraction of H_{2}O that is ice (*s*_{m} + *s*_{ice} = 1) and *X*_{wm} is the fraction of water in mud.

Momentum conservation. Darcy flow in the permeable core bed is as follows(2)where *u*_{m} is the mud velocity and *k* is the permeability.

Substitution of Eq. 2 into Eq. 1 yields the following pressure equation(3)

Permeability. As coarse particles fall from the melting front and settle toward the center, they form the core bed. Until the packing of coarse particles reaches a maximum at a computational mesh grid node, coarse particles may be lofted by upwelling mud convective currents. Once maximum packing is reached at a grid point, we assume that the coarse particles become immobile and form part of the core bed. The remaining mud, consisting of water plus fines, can flow through the core bed. The permeability of each grid cell in the porous, permeable core bed depends on porosity and the size distribution of particles in each grid cell. The effective permeability *k*_{eff} = Π *k*_{i}^{υi}, where *k*_{i} is the permeability of a layer having just the particle size *i* and υ_{i}, is the volume fraction of coarse particles occupied by particle size *i*, and the product is taken over the particle sizes present in each grid cell of the core. There are a number of relationships between porosity and particle size and permeability; here, permeability *k*_{i} is taken from the data of Berg (*41*). Different soil types—sandy, silty, and clayey—are reflected in the particle size distributions for each type. The distribution of this study is a combination of sandy and silty types for coarse particles and silty-clayey for the fines.

Momentum conservation in the mud mantle. Navier-Stokes equation for momentum conservation is as follows(4)

Differentiation of Eq. 4 yields a pressure equation (*42*)(5)

Velocity components are computed from Eq. 4. For these simulations, the terms on the left side of Eq. 4 are very small and are dropped.

Bulk density. (6)where *R*_{x} is the fraction of olivine remaining, that is, the fraction of olivine (and pyroxene) that has not yet converted to serpentine.

Energy transport. Energy transport includes advection of enthalpy plus diffusion plus latent heat conversion plus heat of reaction from serpentinization and radiogenic heating, captured in the energy conservation equation(7)where total energy density is(8)

Thermal conductivity. In the mud and core, *K*_{T} is a volume weighting of the thermal conductivities of water, ice, and grains. In the mantle, thermal conductivity is that of mud times the Nusselt number, as a means of capturing sub-gridscale turbulent mixing. In the frozen crustal layer, *K*_{Tice} is modified by the formulation of Rayleigh number dependence from the study of Mitri and Showman (*43*). Thermal conductivity of ice (*K*_{Tice}) and water (*K*_{Twater}) and specific heat of ice are temperature-dependent and given by Grimm and McSween (*4*)(9)(10)(11)

Thermal conductivity and specific heat of rock are assumed to be constant. Thermodynamic properties of water (density, internal energy, enthalpy, viscosity, and specific heat) are functions of temperature and pressure and are interpolated from tables generated by the SUPCRT92 model (*44*).

Salt transport. (12)In the mantle and ice shell layers, *D*_{c} is enhanced by the same Ra number factor for the thermal diffusivity to reflect turbulent sub-gridscale mixing. A similar transport is applied for *X*_{wm} and *R*_{x}.

Serpentinization. (13)For this reaction, net volume increase in mineral phases approximately equals net volume loss in water, for density of olivine = density of pyroxene = 3200 kg/m^{3} and density of serpentine = 2486 kg/m^{3}, which are the values we used for those minerals. Reaction rate *S* depends on temperature (and presence of water)(14)

Reaction ends in a cell when either water is entirely consumed or olivine and pyroxene are consumed

(15)

Boundary conditions. At the surface, *P* = *C* = 0, *T*_{srf} = *T*_{polar} + Δ*T* cos(θ); θ = 0 at equator. At the surface, *P* = 0, *T* = surface temperature, ∂*C*/∂*r* = 0, and *ū*_{m} = 0 are required. At the mantle-core and mantle-crust interfaces, continuity of *P*, *T*, *C*, and *ū*_{m} is required. At the center, *ū*_{m} = 0, ∂*P*/∂*r* = 0, ∂*T*/∂*r* = 0, and ∂*C*/∂*r* = 0 are required.

Numerical algorithm. The numerical procedure consists, at each time step, of solving for pressure using an implicit numerical discretization; determining velocities, followed by updating mass, energy, and solute transport; and then using the equation of state to reconcile ice and salt (if any) and temperature variables. Time step size is computed after each time step and is controlled by the Courant restriction.

Particle transport. A new feature is the transport of particles. A distribution of coarse particle sizes having different settling velocities can be tracked [although the model allows us to track a range of chondrule/particle sizes, as noted above, if the body accretes a sorted population (the canonical assumption), *T*_{peak} and *T*_{avg} are qualitatively similar, but no further sorting occurs]. For the unsorted simulations, particles are divided into coarse (five sizes from 1 mm down to 10 μm) and fine (<1 μm) classes. The Stokes settling velocity for each particle size is added to the fluid velocity to determine each particle’s total velocity. Fine particles are assumed to remain suspended indefinitely; this assumption is justified for a duration of a few million years of the simulations analyzed. Water plus fines constitutes a mud.

Particle settling rates. Stokes settling relationship governs particles(16)where *V*_{s} is the settling velocity (cm/s), *d* is the particle diameter (mm), ρ_{p} is the particle density, and ρ_{f} is the fluid density.

In our model, particle velocity *V*_{p} = *V*_{s} + *u*_{m}, where *u*_{m} is the mud velocity. At larger sizes (>0.25 mm), particle settling velocity deviates from the Stokes settling relationship. Rubey’s curve provides accuracy in those cases (*19*, *20*).

Settling rates of the sand-sized particles (1 mm down to 50 μm) range from about 2 km/year to about 25 m/year, respectively, in the smallest body considered, which are much faster than the rate at which the thaw front moves outward (~1 m/year). The settling rate for the silty-clayey component (10 μm and smaller) is in the range of 1 to 0.001 m/year, which is on the order of, or less than, the mud convective flow rates, certainly during the most intense phase of convection. The smallest particles can remain suspended for a very long time, especially if flow is occurring. Mud rapidly evolves into a mixture containing mostly the smaller particles.

Porosity. Porosity of grid cells containing particles has a lower limit, specified in the input file. Minimum porosity is also a function of pressure. Porosity ε relates to confining pressure *P* through the relationship(17)where ε_{0} and *p*_{scl} are user-supplied. This relation roughly approximates data on changes in porosity seen in sandstones (*45*).

Viscosity. Viscosity μ_{m} of the mud in our model follows a power law formula(18)where μ_{w} is the water viscosity (a function of temperature), *c* is the volume fraction of particles, and *c*._{lim} is the concentration at which the mixture locks up (*c*._{lim} = 0.75), following Roscoe (*21*). Shoji and Kurita (*22*) used an almost identical relationship. Zoporowski and Miller’s (*23*) viscosities used for mud volcano flows lie in the range of what our model used. Data from the oil and gas exploration literature bracket Roscoe’s curve (*21*). Fritz and Moore (*24*) and Grim (*25*) presented data on sand and clayey particle types. They found roughly one to two orders of magnitude difference between curves for the same concentration. The Roscoe curve (*21*), based on glass spheres, lies roughly midway between the sand and clay curves, about an order of magnitude removed from each. We used it as a compromise between pure sand muds and pure clay/silt muds. In our simulations, the highest initial value of *c* is 0.625. Flow rates and shear rates are very low, and shearing effects are therefore neglected. Viscosity can be affected by other factors not included explicitly in our viscosity model, such as the shape of particles and the particle size distribution; these can also affect the lock-up concentration [for example, see the studies of Delaney *et al*. (*46*) and Stickel and Powell (*47*)]. Once coarse particles reach the core region and the 50% packing density is reached, the core bed structure becomes immobile, and its viscosity is effectively infinite. However, the remaining mud phase (water plus fines) can flow through the core bed’s pores, given sufficient driving force. In the larger bodies, settling rates are faster but convection is also stronger, so the picture remains qualitatively the same.

Core formation and evolution of mud. The core formation process is not instantaneous; it occurs gradually over a period of at least 3 × 10^{4} years and as much as 10^{5} years for some cases. As particles are released at the melting front and fall, their trajectories can be modified by the mud flow occurring around and below them. Mud does not have the same constitution everywhere; in some places, the particulate load is high, and in other areas, it is low, depending on settling history through each grid node.

Convection in the core. When the coarse particles settle into the core and reach their packing density, they are no longer considered part of the mud and cannot move; their effective viscosity is infinite. We then have an essentially two-phase domain. The locked-up coarse particles form a porous, permeable bed that grows in depth until the melt front has reached its limit near the surface. The mud in the pores of the core bed then consists of water plus the fine particles. This mud can still flow through the pores of the core bed, if there is sufficient driving force. The driving force is provided by the energy released during serpentinization and from decay of short-lived radioisotopes. The structure of the core itself, formed by the coarse particles, does not convect (although coarse particles on the topmost layer of the core may still be loftable, until buried by continued raining down of particles from the melt front). Convection through the pores of the core will depend on the permeability of the core bed, in part. Permeability of the core bed is determined from the coarse particle sizes present at each node in the core and the porosity at that point. The mud in the mantle layer consists of water plus fines plus coarse particles that have not yet reached the core region and locked up. Eventually, after the melt process has ended, the mantle mud will also consist essentially of just water plus fines. It can flow given sufficient driving force, but eventually, after several million years for the smaller bodies, the asteroid will freeze completely.

Clustering of particles in the core. There is a burst of convection when serpentinization occurs, very soon after thawing. This energy is being released and drives convection even as particles are flowing down, so there is an opportunity for the flow to change the trajectory of falling coarse particles, especially the smaller ones. Convection will strengthen as the thawed region expands (convection varies with the cube of the domain depth). Upwelling mud currents will entrain and push against the falling coarse particles to varying degrees, imparting to them a horizontal velocity component. As long as the coarse particles are free to move, their trajectories will tend to reflect the convective flow pattern. This results in the clustering of particles seen in the images in Fig. 1. After the serpentinization phase ends (short time frame), and short-lived radionuclides are mostly spent, convection diminishes and the clustering process ends. In Fig. 1, the clustering is most pronounced at slightly more than mid-depth in the core region. This corresponds to the time of maximal convection strength.

Crustal stability. The numerical simulations described in the main text produce density structures that include a frozen crust of variable thickness, generally 2 to 10 km in thickness, during the duration of a few million years of the simulations. The density structure generally shows that the mantle, if it develops, has a lower density than the crust. In that event, questions naturally arise as to the stability of the crust. Formisano *et al*. (*48*) recently analyzed the stability of the crust of Ceres for a variety of conditions. We do not believe that crustal stability is an issue in our study for the following several reasons: (i) Formisano *et al*. (*48*) considered long-term crustal stability for Ceres, after core formation and after short-lived radionuclides have decayed. We are concerned with early-time behavior, in the first few million years, when most particle sorting would likely occur. Even for a Ceres-sized case, they found that crustal stability holds for at least 100 My in all the cases they included. (ii) The mantle density in all cases that we have examined is greater than or equal to the 1440 kg/m^{3} used in the stability analysis of Formisano *et al*. (*48*), meaning that the density difference between crust and mantle would be no more than in their cases. (iii) In the smaller cases we consider, the crustal thickness diminishes for a short period of time to only 1 to 2 km, especially for the non–core-forming cases. During that time interval, flow would mix the near-surface environment. (iv) According to Formisano *et al*. (*48*), crustal stability holds as long as viscosity exceeds a critical viscosity value. For the smaller bodies that we are considering, the stability criterion will be weaker, because critical viscosity of Formisano *et al*. (*48*) varies directly with gravity. For our 50- and 100-km-sized bodies, the critical viscosity is roughly 10 and 5 times lower than that for Ceres, respectively, meaning it is more likely that these smaller bodies will have a stable crust. (v) Desiccation of the surface will reduce the density of the crust, further dampening the Rayleigh-Taylor instability mechanism. In 100 My, for example, a crust several kilometers in thickness could likely lose a significant amount of its H_{2}O ice through sublimation and diffusion through the crust. (vi) In the Rayleigh-Taylor stability analysis, the basis for the Formisano *et al*. (*48*) analysis, perturbations grow at a rate that is proportional to exp(√g) so that, for the smaller bodies, the rate of growth of any crustal perturbation will be very much slower than that for a Ceres-sized body.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

**Acknowledgments:**

**Funding:**This work was funded by the Australian Research Council via their Laureate Fellowship program.

**Author contributions:**P.A.B. and B.J.T. contributed equally to this work.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).