Large Carbonaceous Chondrite Parent Bodies Favored by Abundance–Volatility Modeling: A Possible Chemical Signature of Pebble Accretion

Primitive meteorite groups such as the Vigarano, Mighei, and Karoonda carbonaceous chondrites have enigmatic patterns of elemental abundances, with moderately volatile elements—those that transition from vapor to condensate between ∼400 and ∼900 K—defining plateaus of subequal abundances despite a wide range in volatility. In detail, each group defines a plateau with distinctive nonmonotonic “chemical fingerprints” that have been attributed to combinations of mixing, vaporization/condensation, and fluid-mediated metasomatism—but the extent to which these processes can reproduce the observed variability has not been quantified. Starting with primitive Ivuna chondrite, a two-stage, two-component equilibrium condensation–vaporization model—with gravity implemented as Jeans escape—can explain large-scale plateaus in these chondrite groups, as well as more complex, nonmonotonic small-scale variations. For all three chondritic meteorite groups, models favor earlier high-temperature fractionation under low-gravity conditions followed by a low-temperature fractionation event that took place on a protoplanet at least as large as Ceres. The second fractionation event may represent the fractionation of incoming materials to the planetesimal during protracted pebble accretion. Models with only thermally driven volatile loss, gravity, and mixing can explain more than 80% of the observed compositional variability in these meteorite groups. In our five-parameter model, using only five randomly selected elements yields uselessly large ranges of planet sizes and temperatures, ranges that converge with increasing numbers of elements. These results suggest that even simple models are prone to generating inaccurate conclusions when constrained by too few observations, a fault likely held by more complex models as well.


Background
Carbonaceous chondrite meteorites (CCs) are primitive rocks that represent some of the least processed materials to have survived from the earliest portion of solar system history.The most primitive CCs-the Ivuna or CI chondrites-are thought to be the most geochemically representative of the starting composition of the solar system, based on their similarity to elemental abundances in the solar photosphere (Lodders 2003).The Karoonda (CK), Vigarano (CV), and Mighei (CM) groups have similar but distinctively different chemistries and are thought to have been generated from similar primitive starting materials (Larimer & Anders 1967).
The similarities and differences between the CCs can be demonstrated by comparing abundances (Lodders & Fegley 1998;Lodders 2003) normalized to CI and a refractory element (e.g., Al).Normalized enrichment/depletion profiles of CM, CV, and CK chondrites (Figure 1) all show an overall trend of decreasing abundance toward lower 50% condensation temperatures (T 50% ), as would be expected if more volatile elements are depleted to a greater extent than more refractory elements.However, the decrease in abundance with increasing volatility is not monotonic, as would be expected if initial abundance and volatility were the only factors.It has been noted previously that the moderately volatile elements (MVEs) with T 50% between ∼400 and ∼900 K form a "plateau" where the overall slope is near zero.This plateau has been previously explained as the product of two-component mixing between a primitive component and a component that is depleted in volatile elements (Larimer & Anders 1967;Palme et al. 2014;Braukmüller et al. 2018).
Although the MVE plateau can generally be described as level (i.e., having a slope near zero), it is not flat: local variations in normalized abundances between neighboring MVEs with similar T 50% exceed a factor of two.These local variations or "chemical fingerprints" are unique to each meteorite group and are somewhat enigmatic, as they are not easily explained by volatility-controlled processes such as vaporization or condensation.
Vaporization of an element does not necessarily equate to loss from a body to space, and not every element or molecule of vapor faces the same hurdles in their attempted escapes (Lammer et al. 2020).In differentiated planetary bodies, the formation of a core results in siderophile and chalcophile elements being concentrated near the center of the body, where escape from the surface is difficult at best, whereas the formation of a crust would concentrate other elements close to the surface, making escape as vapor more feasible.Even in the simplified case of a homogeneous, well-mixed spherical body, once each element is at the surface and in the vapor phase, it must be able to overcome the gravitational attraction of the object from which it arises.One simple parameterization of the effect of gravity is Jeans escape (Jeans & Jeans 1952;Young et al. 2019;Tang & Young 2021): particles with velocities (v i ) greater than the escape velocity (v e ) of the body have a higher likelihood of escaping than particles with v i < v e .At a given temperature, the fraction of vapor particles with v i v e is lower for more massive vapor species than for less massive species.The mass and T 50% of an element are uncorrelated, and as such Jeans escape is a plausible mechanism for fractionating elements with similar T 50% provided that they have sufficiently different masses.The least complex models of vaporcondensate fractionation assume equilibrium between the condensed phase (the object or body) and the vapor phase (the atmosphere or nebula).In single-stage equilibrium scenarios, loss by vaporization and gain by condensation are indistinguishable, and Jeans escape can result in preferential loss of less massive species or preferential retention of more massive species.

Vaporization-Condensation Model
To determine the extent to which a simple model (Figure 2) of equilibrium vapor-condensate fractionation combined with mixing can explain observed chemical fingerprints, we developed a numerical model of the formation of CCs, where element volatility is calculated by using the Clausius-Clapeyron equation to extend published T 50% values (Lodders 2003;Sossi & Fegley 2018;Wood et al. 2019) to arbitrary fractions of volatile loss.Unlike previous efforts to explain CC sample chemistry through multicomponent mixing of known materials, we begin with only primitive CI chondrite (horizontal line at y = 1 in Figure 1) and generate the components using temperature-dependent volatile modeling, rather than using prescribed reservoirs.
Starting with CI chondrite, the equilibrium vapor and condensate fractions are calculated for a fractionation event with a characteristic T and v e .For each fractionation event, only the fraction of the vapor with v i v e is "lost."The portion retained is the fraction of the vapor with v i < v e and the condensed fraction.This retained portion is then combined with CI in varying proportions.The mixture is subjected to a second fractionation event (FE B ) at any of a range of T and v e , and again the condensed plus low-velocity vapor portion of the mixture is retained.This final composition is then compared to observed CC compositions, with the agreement quantified by the adjusted R 2 between model output and measured values, all normalized to CI chondrite and the abundant refractory element Al.Using a forward-modeling, model-space search approach, more than 10 6 combinations of parameters (T A , T B , R A , R B , mixing ratio) are compared to the observed element abundances of each CC to map out best-fit regions of parameter space.

Results
Fractionation models only considering the MVEs with T 50% between 400 and 900 K (N = 18-20 elements) generate relative abundance-T 50% relationships that are remarkably similar to the complex patterns of the chemical fingerprints observed in CC chemical compositions (Figure 3).Based on the adjusted R 2 between the observed and model MVE abundances, the best-fit conditions from our model explain 81% of the variation in MVE abundances in the CVs, 50% of MVE variations for the CMs, and 70% of the variation in CK MVEs.
Qualitatively, our model reproduces the distinctively flatbottomed "bucket" shape formed by the sequence Rb-B-Zn-Se-Ag, with subequal abundances of Zn and B (Figure 3).Models also reproduce the relative abundance fingerprints of Cd-Pb-In-Bi-Cl-Br in CMs and Pb-In-Bi-Cl-Br in the CVs.Ratios of refractory elements over their more volatile neighbors, such as K/Au, S/Te, and Sn/Cs, are also less than unity for all three groups, as observed.Models of CVs, CMs, and CKs can replicate the observed patterns for K-Au-Sb-Ge-Rb despite the fact that the patterns are not all the same in the three CCs studied here (i.e., Rb > Ge for the CVs and Ge > Rb for the CMs and CKs).
Fitting models with additional elements that are in between MVE and refractory elements (e.g., up to As with T 50% = 1235 K) results in statistically improved fits (R 2 = 0.62-0.86)with similar qualitative fitting characteristics.Excluding elements more refractory than Ba (T 50% > 1423 K)-which are not the subject of this study-these models explain 81% (CV and CM) to 84% (CK) of the variation observed.
Using a model-space search approach (Beghein 2010; Xing & Beghein 2015), we can explore large regions of parameter space (temperature, planetary body size, and mixing ratio) to find regions of model space where the data are best fit (Figure 4).In general, best-fit mixing ratios of high-and low-T components (CK = 1.9, CM = 0.6, CV = 1.5) are consistent with the nonmatrix (high-T) to matrix (low-T) ratios derived by previous workers from petrographic or geochemical data (McSween 1979;Hewins et al. 2014;Braukmüller et al. 2018;Alexander 2019).
Best-fit models also consistently have the earlier, high-T fractionation event on smaller bodies/objects and the later, low-T volatile loss event on a larger body (Figure 2).CV, CM, and CK chondrites are best fit with nearly identical early heating conditions of T > 1000 K, as objects no larger than the size of the asteroid Vesta.This is an upper limit in object size and is consistent with the observed dimensions of various CC components that underwent high-T geochemical processing  (Sossi & Fegley 2018), normalized to CI and to Al and plotted against T 50% from Wood et al. (2019).Note that although all three meteorites plotted here have a low-slope MVE "plateau" between 400 and 900 K, there are local variations between neighboring elements with similar T 50% but abundances that differ by up to a factor of two.
Constraints on the second fractionation event are also similar for each of the three meteorite groups: for the CVs, the data are best fit by models with T = 800-950 K on a body larger than Ceres (radius > 500 km).Data for CM chondrites are best fit by the same body size as for the CVs but slightly lower temperatures (700-900 K).Modeling the CKs also yields similar results, with the best-fitting models having both size (larger than Ceres) and temperature (750-950 K) constraints overlapping the CV and CM groups.
For all three CC groups studied here, the second fractionation events are modeled to have occurred at temperatures similar to the highest temperatures indicated by petrographic or geochemical thermometry of the chondrites themselves, though this agreement is abetted by the large range of temperatures reported for each of the different CCs: for the CM chondrites, peak metamorphic temperatures range from <600 to over 1000 K (King et al. 2021;Velbel & Zolensky 2021).Temperature estimates for the CV chondrites also vary considerably: various authors (e.g., Brearley & Krot 2013;Bonal et al. 2020) have noted the wide range of temperatures (330°C-700°C; 603-973 K) for the CV Allende alone.Our best-fit temperatures for CVs (800-950 K) are consistent with the upper end of the temperature range from the literature.Peak metamorphic temperatures for the CK chondrites vary a great deal depending on which CKs (e.g., Righter & Neff 2007;Chaumard & Devouard 2016) are selected, but the range is generally consistent with our model results (750-950 K).The agreement between our model-derived best-fit temperatures for the second fractionation event (FE B ) and previously determined constraints on peak parent-body metamorphic temperatures could be taken as evidence that they are one in the same.However, this is not required in our model.
In addition to thermal losses, the aqueous alteration experienced by all of these meteorites may or may not have been isochemical (Bunch & Chang 1980;Brearley & Krot 2013;Piralla et al. 2021).The goal of this study is not to dispute that hydrothermal metasomatism occurred, but rather to ask what fraction of the observed chemical fingerprints can be generated with mixing between components generated by equilibrium fractionation between refractory and volatile components.This is an important question to raise, especially given that most CCs show no evidence of fractionation between water-soluble and water-insoluble elements (Bland et al. 2009).For the CV and CK chondrites, 70%-80% of MVE variations can be explained without hydrothermal metasomatism.For the CM chondrites, our model can explain only 50% of the MVE variance.Although the R 2 value for CM models rises to >0.8% if a greater number of elements are considered (T 50% between 400 and 1423 K), it is also possible that the inability of the model to better replicate the MVE chemical fingerprints of CMs is the result of the effect of aqueous metasomatism.It is reasonable to suggest that the observed variations are the product of a combination of mixing, volatile loss, and aqueous metasomatism.As an end-member, our model of equilibrium vapor-condensate fractionation combined with Jeans escape quantitatively explains the majority of the MVE variations observed in each of the CCs studied: if quantitative aqueous geochemical modeling could be performed in a similar manner, it would likewise be an important contribution to our understanding of how the CCs evolved.

Planetary Body Size
Our modeling of the chemical variations observed in CM, CV, and CK chondrites is most consistent with the latter stages of their respective chemical evolutionary histories taking place on bodies that were at least as massive as Ceres, if not larger, as models with smaller second-stage bodies do not generate the fractionation patterns measured in these CCs.As parent-body size increases, the interiors of these bodies become sufficiently hot to undergo differentiation, in apparent contradiction to the canonical model of chondrite parent bodies based on undifferentiated meteorites.This conflict can be resolved by noting that it is possible that at least some chondrites are samplings of undifferentiated material from closer to the surface of a parent body that has experienced melting and differentiation at depth, similar to the model proposed by Elkins-Tanton et al. (2011), on the basis of the magnetic properties of some chondrite samples.An outstanding complication of this model is that we have not yet identified all the different lithologies that would have existed in the differentiated portions of CM, CV, and CK parent bodies-or we have sampled them, but we have not yet come to an agreement on how the pieces of this particular puzzle fit together.However, not all the components of the CCs were entirely created at low pressures: Hamilton et al. (2021) used petrologic evidence to constrain the size of the parent bodies of two primitive meteorites: the CV chondrite Allende, and an ungrouped CC clast in the Almahatta Sitta ureilite (known as Primitive CI chondrite is subjected to elemental fractionation during an earlier fractionation event (FE A ), and the more refractory residue is mixed with additional, unfractionated CI chondrite.This mixture is then subjected to a second fractionation event (FE B ), and the final composition is compared to observed abundances.fragment 202).In both cases, the mineralogy and petrography of the CCs require their respective (and distinct) parent bodies to have been at least 640-1800 km in diameter, consistent with the parent-body sizes inferred here from the MVE compositions of CV, CM, and CK parent bodies.
The influence of the mass of a terrestrial planetary body on its chemical evolution has been highlighted by other recent interpretations of both isotopic and abundance data from a variety of planetary materials: Day & Moynier (2014) attribute the greater volatile content of Earth and Mars-as compared to smaller bodies such as the Moon and Vesta-to their much larger size.Day et al. (2020) suggest that both Zn abundances and δ 66 Zn may be in part controlled by the size of the body, with smaller bodies experiencing greater loss and isotopic fractionation.Tian et al. (2021) note correlations between surface gravity and geochemical parameters such as δ 41 K and Mn/Na and attribute these relationships to planet size.Several groups have pointed out correlations between bulk silicate H 2 O abundances and parent-body radius of rocky inner solar system bodies (McCubbin & Barnes 2019; Sarafian et al. 2019;Tian et al. 2021).Jeans escape is a plausible mechanism by which the mass of a planet can influence how it evolves chemically or isotopically.Jeans escape may not be the primary mechanism of mass loss from most protoplanets (Volkov et al. 2011), but the most parsimonious explanation of the agreement between our model and observed chemical profiles is that the influence of gravity is adequately modeled as Jeans escape, and that the size of the body plays an important role in the chemical evolution of CCs-especially the fractionation of elements with similar T 50% .Future work building on these results might include more nuanced treatments of escape (e.g., inclusion of hydrodynamic escape; Hunten et al. 1987;Young et al. 2019;Tang & Young 2021); the effect of heterogeneous solid bodies, such as that due to differentiation; more rigorous thermodynamic calculations for volatility estimates with a dependence on composition and pressure in addition to temperature; and the consideration of characteristic temperature ranges for each component, as opposed to specific temperatures.

Chemical Topographic Evidence for Pebble Accretion?
Our simplified model posits fractionation at the surface of a planet without the complications of mixing, transport from the interior, or other processes that would be required if the body is constructed prior to the loss of volatiles.However, the possibility that planetary bodies can grow to the size of Ceres or beyond by pebble accretion (Lambrechts & Johansen 2012;Bitsch et al. 2015) provides a model by which the CV, CM, and CK parent bodies could have acquired their distinctive gravity plus volatility chemical topography during construction.Incoming pebbles that were partially or wholly ablated by wind erosion (Demirci & Wurm 2020), vaporized during their collision with the protoplanet, or vaporized by the impact of subsequent impactors would deliver material to the parent body's gravitational well at or near the planetary surface.The combined effects of volatility and Jeans escape would then be poised to influence the separation of low-mass volatile elements from elements with higher mass and or lower volatility.This fractionation during pebble accretion hypothesis is consistent with our simple model, which explains a large portion of the fractionation patterns observed in the meteorite groups studied here.If chemical topography is related to pebble accretion, then it should be possible to construct testable hypotheses relating chemical or isotopic measurements of meteorites to pebble accretion conditions, provided that the chemical or isotopic signature of fractionation from accretion can be deconvolved from the effects of multicomponent mixing and other nonisochemical processes.Au) and depletions (e.g., K) in elements with higher T 50% .

How Many Elements Are Enough?
Our five-parameter model should be perfectly constrained by only five elements, but models of natural compositions run with randomly selected subsamples of elements with T 50% < 1000 K indicate that models with five elements can yield best-fit models that vary by more than five orders of magnitude in R.A. and temperatures for each component that span >800 K.At least nine MVEs are required to have a 95% probability of mutually resolving the two fractionation events in temperature and radius using a single set of nine elements (Figure 5).But this does not guarantee that any one set of elements will lead to success: for both natural and synthetic chemical profiles where the temperatures and radii are known or specified, models with too few parameters can be accurate as an ensemble, but with very large variance, meaning that any one individual model fit to the subset of elements selected is likely to yield conclusions that are different from models fit to a different subset of elements.
More complex models with greater numbers of free parameters (e.g., more than two components, more complex and variable loss functions, internal heterogeneity, pressuredependent volatility, complex geodynamic or atmospheric models) will require a greater number of independent constraints (with each element's abundance or isotope ratio serving as a single constraint).Models of planetary-scale chemistry that rely on few constraints (e.g., a single isotope ratio) are at risk of being mathematically underconstrained and unlikely to yield robust results.Future modeling efforts would benefit from testing for sufficient constraint by comparing results of models constrained by the maximum number of observations to those of models constrained by randomly subsampled combinations of elements where N subsample < N.
However, no published set of T 50% values eliminates unexplained chemical topography, so we defer to the more recent values of Wood et al. (2019) for our modeling and reserve the values of Lodders (2003) as an independent test of our model as described below.The complex pattern of MVE variation between the three CCs depicted in Figure 1 is similar but not identical.This similarity argues that these patterns are robust features of the CCs and not the product of random error in the measured elemental abundances of chondrite meteorites.The fact that the topographies are not identical-a statement that is true regardless of which set of T 50% values are chosenargues against this simply being the result of errors in the calculated T 50% values.Models of MVEs where the respective T 50% values have been randomized have no apparent information content and do not yield results that are consistent with models that use established T 50% values, or even with other randomized trials.

Appendix B Model Overview
Designed as a test of the two-step model of CC parent-body formation, the model presented here (Figure 2) starts with CI composition material and calculates the residue present after loss at some specific temperature (T A ) and planet size (R A ).This material is then combined with unheated CI composition material (at varying mixing ratios m).
The residue of this mixture is then calculated for loss at a second set of temperatures (T A ) and planet sizes (R B ).For wide ranges of T A , R A , T B , R b , and m parameter values, the resulting compositions were compared to geochemical observations, with regions of minimum misfit representing the best-fit portions of parameter space.A model-space, forward-modeling search approach was used to find the best-fit model by misfit minimization, for a wide range of values for each parameter (with a total of >10 6 parameter combinations tested for each CC group).

Appendix C Volatility Model
The 50% condensation temperature (T 50% ) is a metric of volatility that represents the temperature at which there are an equal number of moles in the vapor and in a condensed phase for any given element, under specific conditions such as those relevant to the solar nebula (Lodders 2003).Beginning with the T 50% values of Wood et al. (2019), we extend T 50% to arbitrary vapor fractions as a function of temperature (T X% ) by asserting that the vapor fraction would scale with vapor pressure, which in turn can be related to temperature using an equation derived from the approximate form of the Clausius-Clapeyron equation: where P is the vapor pressure, T is the temperature, L is the specific latent heat, R is the specific gas constant, the subscript X represents an arbitrary vapor mole fraction, and the subscript 50 indicates the mole fraction of vapor at the reference T 50% from Wood et al. (2019).Rearranging so that which is a linear equation relating ln(P X ) to 1/T X , for any one element.If we assume that the slope and intercept (square brackets in Equation (C3)) are both constants for any given element, we still need one additional constraint because one cannot solve a linear equation with only one parameter (T 50% ).
We set up a second set of equations similar to Equation (C3), but with T 100% instead of T 50% : The temperature at which 100% of the moles of a given element are condensed (T 100% ) is absolute zero for every element.But as dividing by zero is not preferred, we instead define T 100% * as the temperature at which the vapor pressure is effectively zero.For any given criterion of "effectively zero vapor," it is reasonable to assert that T 100% (and therefore T 100% * ) would be related to T 50% in the sense that elements with high T 50% should have higher T 100% than would elements with low T 50% .We use a single tuning parameter (tau) to relate T 100% * to T 50% for every element, regardless of the number of elements selected: We can then solve a series of linear equations for the slope and intercept parameters in Equations (C3) and (C4) for each element.It is important to note that the vaporization model only has a single value of tau for all the elements and temperatures in any given model, regardless of the number of elements or temperatures modeled.For N elements and an arbitrary number of temperatures, the model requires N + 1 constraints: T 50% for each element, and tau.
In order to test this simplified evaporation model, we can use it to predict the T 0% (100% vapor) and then compare model values to those calculated by Lodders (2003).For t = 33, the model accurately fits all the elements for which both T 50% and T 0% were provided by Lodders (2003).The mean deviation between the lines best fit to the Lodders (2003) data and our model is 4.8 K, decreasing from ∼7 K below 90 K to <3 K above 1550 K.This represents 1.5% error for elements with T 50% > 400 K.This same value of t also yields model T 0% that matches the T 0% of H 2 O derived from the models of Wagner et al. (2011).
The tuning parameter tau controls the curvature near T 50% : large values of tau result in a more gradual increase in vapor fraction, whereas smaller values of tau make the transition from ∼zero vapor to 100% vapor closer to a step function.Changing the steepness of the curves results in only a very small effect in the final model fitting.Regardless of the temperature range of the transition from mostly condensed to mostly vapor, the model retrieves a very similar result, with the major difference being that the best-fit regions for step function models are smaller than for thermal models with larger tau.This suggests that geochemical models with two components and two thermal events can use a step function vapor fraction versus temperature model to derive useful results.This is not to say that experimental constraints or calculations are not important, just that the overall behavior of each element is sufficiently constrained by T 50% .
There are multiple formulations of T 50% in the literature (Lodders 2003;Fegley & Schaefer 2010;Wood et al. 2019), and although these formulations are highly correlated, there are elements for which one estimate is significantly different from the others.As a test of the extent to which changing T 50% changes the results, models based on Wood et al. (2019) were compared to models based on Lodders (2003).For CV chondrites, best-fit models with evaporation defined by these two sets of T 50% are very similar.This is because the models are constrained by many elements, most of which are not the subject of much disagreement between the various authors.
Our model allows us to calculate the vapor fraction of any element at any temperature, provided that the remainder of the conditions are consistent with those for which the original T 50% values were derived.In theory, this approach can be used to extrapolate the fraction of condensate using any data set (e.g., Lodders 2003;Fegley & Schaefer 2010) and starting from any percent (e.g., the T 99% of Tartèse et al. 2021), provided that the correct tau value can be determined.
The accuracy of our volatilization model hinges almost entirely on T 50% values, which in turn are specific to the conditions of equilibrium between condensed and vapor phases.However, our model best-fit temperature ranges correspond to nebular conditions with low pressure and low oxygen fugacities.Sossi et al. (2019) demonstrate that increasing pressure and oxygen fugacity decreases the vapor fraction at a given temperature, resulting in higher T 50% .If the conditions during the evolution of these CCs were more oxidized (or at higher pressure, or both), the model would yield even higher temperatures for both events.For example, (unwisely) using the values of Sossi et al. (2019) in our model would yield best-fit temperatures for CCs that are more than 800 K higher than temperatures derived from models based on the preferred T 50% of Wood et al. (2019).

Appendix D Jeans Escape
The approach described above provides a model of vapor fraction at any temperature for any element.However, as vaporization does not necessarily equate to loss, we add a Jeans escape parameter to the evaporation model to take into account the gravitational force opposing the loss.Here we frame the discussion by interpreting the gravitational potential in terms of the size of a protoplanetary body from which a given element is attempting to escape, but other scenarios can be envisioned (e.g., escape from the proto-Sun's gravity during the disk phase of solar system evolution).The escape velocity of a planet is given by where G is the gravitational constant, M is the mass, and r is the distance from the center of mass (here assumed to be the radius of the planet).
At any given temperature, the distribution of particle velocities in a vapor can be modeled as a function of temperature and vapor species mass by the Maxwell-Boltzmann relationship.For a given reference velocity (here v esc ), the probability P v of a particle having a velocity greater than that of the reference velocity is given by where a is defined as with m the mass of the particle, T the absolute temperature, and k B the Boltzmann constant.As one would expect, higher temperatures and smaller particle masses yield distributions of velocities that are greater than at lower temperatures or for more massive vapor species.For a given vapor species and temperature, a fraction of the particles will have velocities greater than the escape velocity of the body.More massive planets have higher escape velocities and retain a greater fraction of any volatile species.The choice of Jeans escape over other escape mechanisms that might be operating is a pragmatic one, as our goal is to add realistic mass-dependent chemical topography to the model with as few parameters as possible.
This results in what is effectively a model of volatile loss from an isothermal, well-mixed spherical body as a function of mass and temperature.It should be noted that the entire body need not be well mixed, only the portion of the body for which the observed compositions are representative.In order to account for the existence of other mechanisms of volatile loss that are not sensitive to gravity, we decreased the magnitude of the Jeans escape filter to 50%.This tends to decrease the magnitude of the local chemical topography in the model output-though in many cases the amplitude of the modeled local topography still exceeds the observed topography.It should be noted that all model tests with Jeans escape >0% fit the data better than models that have 0% Jeans escape but are otherwise identical.

Appendix E Vapor Speciation
The choice of vapor species changes the behavior of elements in the model because the velocity at a given temperature is a function of the mass of that vapor species.Consider as an example a metal (M) where the most likely species are the metal M x+ or the oxide of that metal (e.g., MO x/ 2 or M 2 O x ).Because the oxide is more massive, a smaller fraction of the distribution of velocities at a given temperature will be above the escape velocity than would be for the metal ion without oxygen.Different elements make different oxides -and other vapor species-and as such changing the vapor species can change elemental fractionation, even between elements with similar T 50% .As an example, the elements Zn and B define a small range of T 50% (∼700-740 K) and a large range of masses (as metal ions) from ∼11 to ∼65 AMU.As metals, the loss of Zn would be favored over B because of its lower T 50% , but the higher mass of Zn makes it less likely to overcome the escape velocity of the planet at a given temperature.However, as an oxide, B could form B 2 O 3 , which has a molecular mass of ∼70 AMU, a much greater change than for Zn, which only sees its mass increase from ∼65 to ∼81 AMU when oxidized to ZnO.We might therefore expect greater depletions of B (relative to Zn) if the vapor species were B and Zn, as opposed to B 2 O 3 and ZnO.Models of volatile loss run with the fraction of oxidized species for each element varying between 0% and 100% follow expected patterns, where greater proportions of elemental B result in greater depletions of B relative to Zn.The model is capable of generating B > Zn and B < Zn relationships as a function of oxidation state, consistent with observations of both B > Zn and B < Zn in CCs.It may be that this pair of elements (and other similar sets with similar mass and T 50% relationships) could be used to evaluate the relative importance of oxidized species.However, given that models with two-component mixtures are capable of generating B > Zn or B < Zn without changing volatile species oxidation state, B/Zn ratios of complex natural samples may be difficult to interpret in terms of oxygen fugacity.

Appendix F Refractory Elements
Our model is neither intended to replicate nor capable of replicating the striking topography observed in highly refractory elements (T 50% > 1400 K; Figure 1).That a simplified, two-component model cannot completely account for the compositions of the CCs studied here is unsurprising and in agreement with previous work.Alexander (2019) noted that refractory inclusions in CCs have enrichments (e.g., W vs. Re and Os) that are consistent with observations in the CVs, CKs, and CMs.However, other enrichments, such as Zr/Hf, are not consistent between these meteorite groups, with Zr/ Hf > 1 in the CMs and CVs but Zr/Hf < 1 in the CKs.This is difficult to explain with the addition of any third component common to all CCs.Even the fourth H 2 O-rich component postulated by Alexander (2019) is an unlikely carrier of either Hf or Zr.
Our model is a simple but effective framework that satisfactorily emulates some models of planet formation (i.e., pebble accretion) and can be further developed to better account for additional processes and components.This development might be in the physics of the model, the limit of two components, or allowing for models that incorporate recondensed vapors-not only refractory residues-as potential components.

Appendix G Mixing Ratio
Our model also determines the mixing ratio of high-to lowtemperature components that generates a best fit for each CC group.This is one of the most robust constraints provided by our model and is almost entirely controlled by the normalized abundances of the MVEs (Braukmüller et al. 2018).If we assume that our high-T and low-T components correspond to the nonmatrix and matrix components proposed by Braukmüller et al. (2018), we can compare our mixing ratios to previous estimates.Best-fitting nonmatrix/matrix ratios derived from our model (CK = 1.9, CM = 0.6, CV = 1.5) match closely the ratios derived by previous workers from petrography, or modeling of geochemical or isotopic data (McSween 1979;Hewins et al. 2014;Braukmüller et al. 2018;Alexander 2019).This is consistent with-but does not conclusively prove-the idea that the low-and high-temperature components in our model correspond to matrix and nonmatrix components, respectively.

Appendix H Model Testing and Quality Assessment
The ability of the model to retrieve specified temperatures and radii was tested for a wide range of parameter space, with varying amounts of Gaussian noise added.The goodness of fit is demonstrated by the reduced R 2 value, which is a useful metric approximating the extent to which a given model explains the observed variations.
Because this is a five-parameter model, model-space searches require considerable time, especially as resolution increases.Initially, brute-force model-space searches were run in stages with increasing resolution (e.g., from ΔT = 20 K; 10 K; 5 K; 2 K; 1 K) until results were consistent.However, models with ΔT < 10 K were not significantly better fitting or more informative than models run with ΔT = 10 K, and as such 10 K was chosen as the temperature resolution of the model.Planetary body size and mixing ratio resolutions were chosen in a similar manner.

Appendix I Degrees of Freedom and Uniqueness
There are five free parameters in the model, so one would expect that at least five elements would be required for a robust result.Models attempting to fit synthetic chemical fingerprints created for a known pair of temperatures with 1% Gaussian noise were used to assess the extent to which the model behaves according to expectations.Synthetic chemical profiles generated under known conditions require at least seven independent observations to resolve the best-fit values given their uncertainties.However, making seven observations does not result in an accurate model result every time, as individual models with different elements can yield very different outcomes.Synthetic models constrained by only three randomly selected MVE elements (26 elements taken three at a time) always result in very precise fits to chemical profiles, but the best-fit models correspond to within ±50 K of both prescribed temperatures in only ∼8% of cases, 28% for N = 5, and 78% for N = 10.A ∼95% success rate is not achieved until there are at least 16 elements, a requirement that increases to 20 elements for models with 10% Gaussian noise.For input temperatures of 1100 and 700 K, models with 10% noise but constrained by N 20 observations yield an ensemble of bestfit temperatures of 1087 ± 67 K and 700 ± 3 K (2σ), respectively, with a maximum error of only 70 K.If we only use half as many observations (N = 10), then the maximum error rises to 500 K, demonstrating that the elements that happen to be available are critical for small numbers of elements, but for large numbers of elements the absence of random elements is less critical.
Similar observations can be made using subsets of the measurements used in the modeling for this work.As an example, CV chondrite with N = 21 elements (20 MVEs plus Al as a reference element) yields best-fit temperatures of 1500 and 900 K and best-fit parent-body sizes of ∼300 and 2000 km for the early and late fractionation events, respectively.
Omitting one element at random yields identical results, but as the number of elements randomly omitted increases, so does the variation in the best-fit values.Randomly omitting 7 of the 21 elements (13 MVEs + Al) results in temperatures in the range of 860-1590 K (with a mean of 1320 K) and 760-1570 K (with a mean of 909 K), respectively.These averages are still within error of the answers derived from N 20 data sets, but any given set of 13 MVEs may yield best-fit temperatures that vary wildly.Parent-body size estimates are also affected to a lesser extent, with best-fit values in the range of 0.01-300 km (with a mean of 100 km) and 1000-3000 km (with a mean of 2000 km).Various random samplings of elements for CV chondrites are shown in Figure 5, and taking into account the correlations between temperature and size, we see that data sets with N 9 (8 MVEs + Al) are resolved from each other at 2σ.To have 95% confidence of obtaining a value within 100 K of the value derived from N = 21 elements requires N 20 for the high-temperature fractionation, but for the second, low-T event, it requires only N 15 elements.This highlights the lower sensitivity of the model to the early fractionation event that is overprinted and obscured by the later fractionation event.
The results presented here suggest that caution should be applied when interpreting models of planetary geochemistry constrained by a small number of independent constraints because even in those cases where highly precise and accurate fits to the data are obtained it is possible or even probable that the model has failed to retrieve the "true" value.Models with fewer parameters than the model described here may require fewer constraints, but it is difficult to envision a scenario where one or two constraints yields an accurate best fit every time.To the extent that the modeling utilized here to study fractionation of elements applies also to isotope ratios, our work would suggest that even a very precisely determined isotope ratio of a single element does not provide sufficient constraint when simulating scenarios of volatile loss with multiple steps and multiple components.

Figure 1 .
Figure1.Abundances of elements in CV (cyan), CM (brown), and CK (magenta) meteorite groups(Sossi & Fegley 2018), normalized to CI and to Al and plotted against T 50%from Wood et al. (2019).Note that although all three meteorites plotted here have a low-slope MVE "plateau" between 400 and 900 K, there are local variations between neighboring elements with similar T 50% but abundances that differ by up to a factor of two.

Figure 2 .
Figure2.Conceptual overview of the numerical model utilized in this study.Primitive CI chondrite is subjected to elemental fractionation during an earlier fractionation event (FE A ), and the more refractory residue is mixed with additional, unfractionated CI chondrite.This mixture is then subjected to a second fractionation event (FE B ), and the final composition is compared to observed abundances.

Figure 3 .
Figure3.Model fits to abundances of elements in CV (cyan), CM (brown), and CK (magenta) meteorite groups for elements between 400 and 1000 K. Observations (dashed) are compared to the best-fit model (solid black line).The gray region indicates the ensemble of best-fitting models (N = 5000 of more than 5 × 10 6 combinations).The MVE plateaus and the local topographies within are replicated by the model, as are some systematic enrichments (e.g., Au) and depletions (e.g., K) in elements with higher T 50% .

Figure 5 .
Figure 5. Error ellipses for random subsampling of Al plus 5-20 MVE elements for CV chondrite, where [T A , R A ] pairs are plotted in magenta and [T B , R B ] pairs are plotted in blue.Each pair of ellipses represents the 95% confidence regions for [T A , R A ] and [T B , R B ] pairs, for 100 random sets of N = 5-20 elements.Mean values(cyan and black points) show an order of magnitude of scatter in R.A. and 300 K in T A , with a systematic shift toward indistinguishable T A and T B conditions.For a model constrained by only five elements-depending entirely on the choice of elements-each best-fit temperature can vary by a factor of two.However, larger numbers of elements deliver consistent results, and models having N 9 elements are distinguishable in T-r model space at the 95% confidence level.