Thermal Extraction of Volatiles from Lunar and Asteroid Regolith in Axisymmetric Crank–Nicolson Modeling

There is growing interest in extracting water from the Moon (Casanova et al. 2017), from Mars (Abbud-Madrid et al. 2016), and from asteroids (Nomura et al. 2017). NASA’s Lunar CRater Observation and Sensing Satellite (LCROSS) demonstrated the existence of water ice on the Moon when it impacted into the Cabeus crater, a permanently shadowed region (PSR) near the Moon’s south pole. The resulting ejecta blanket was found to contain water and other volatiles (Colaprete et al. 2010; Gladstone et al. 2010). Carbonaceous asteroids contain hydrated and hydroxylated phyllosilicates (Jewitt et al. 2007). These volatiles are stable at typical asteroid temperatures in near-Earth orbits in vacuum, but the water evolves when heated to moderate temperature (Zacny et al. 2016). Mars has abundant water in the form of glacial deposits, hydrated minerals including polyhydrated sulfate minerals and hydrated phyllosilicates, and water adsorbed globally at a low weight percent onto the surfaces of regolith grains (Abbud-Madrid et al. 2016).

The extraction of water on the Moon or Mars can be accomplished through strip mining with subsequent processing of the mined material or through in situ thermal techniques: injecting heat into the subsurface, providing a means for vapor or liquid to reach the surface, and collecting it in tanks. Methods on asteroids could be similar or could involve bagging the entire asteroid and heating or spallation of the rocky material with concentrated sunlight (Sercel et al. 2016). Water can be used to make rocket propellant to reduce the cost of operating in space (Sanders et al. 2008; Hubbard et al. 2013; Sowers 2016; Kutter and Sowers 2016), can serve as passive radiation shielding for astronauts (Parker 2006), and can provide life support (Kelsey et al. 2013). In addition to supporting national space agency activities, water-derived rocket propellant can be used commercially for boosting telecommunication satellites from low-Earth orbit or from geosynchronous transfer orbit into geostationary orbit, or for supporting space tourism or other nongovernmental activities (Metzger 2016).

It is difficult to test extraterrestrial water extraction technologies on Earth because of the high preparation costs for realistic test environments: large-scale beds of frozen regolith in vacuum or Mars atmosphere chambers (Kleinhenz and Linne 2013). Nevertheless, testing is vital, as shown by Zacny et al. (2016) performing thermal extraction of water from icy lunar simulant. Those tests found that water vapor moves through the regolith down the thermal gradient away from the hot mining device, so depending on the particular geometry of the system, it could either collect large amounts of water or none at all. This is highly dependent on environmental conditions including vacuum, ice characteristics, and thermal state, so testing without the correct environment would have little value. However, the environments can be simulated numerically to perform low-cost digital design evaluation in lieu of some of the testing, reducing the cost and speeding the schedule of hardware development. For this approach to work, the equations must accurately predict the thermodynamic behaviors of icy, extraterrestrial soil or hydrated minerals in hard vacuum or in low-pressure conditions and at extreme temperatures. This is still challenging because most measurements of thermal conductivity for soils have been performed in terrestrial conditions with liquid water content, Earth’s atmospheric pressure, and/or ambient temperatures. Therefore, more work is needed to develop improved models.

The application that led to the present effort is the World Is Not Enough (WINE) spacecraft concept, which is being developed by Honeybee Robotics under a NASA contract (Zacny et al. 2016). WINE will be a small spacecraft, approximately 27U in CubeSat dimensions (3 by 3 by 3 cubes), with legs for walking short distances and a steam propulsion system for hopping multiple kilometers (Metzger et al. 2016). WINE spacecraft could operate on a body such as dwarf planet Ceres obtaining water from hydrated minerals that may exist on its surface, or on a moon such as Europa where ice is abundant. WINE will drive a corer into the regolith to extract water and perform science and prospecting measurements on the regolith. The water will be extracted thermally by heating the material in the corer. Vapor will travel into a collection chamber where it is frozen onto a cold finger. After multiple coring operations have collected enough water, the tank will be heated to high pressure and vented through a nozzle to produce hopping thrust. The development and validation of the coring and water extraction system require at least two-dimensional (2D) (axisymmetric) computer modeling of heat transfer in the regolith. The modeling is needed to determine energy requirements for this process to set requirements for the spacecraft power system and to determine whether solar energy is adequate or whether radioisotopic heater units (RHUs) are needed to generate adequate thermal energy on a particular planet.

The authors were team members for another application of this modeling: NASA’s Resource Prospector mission (which was canceled while in development). It was planned to prospect for water in the Moon’s polar regions by drilling into the regolith, bringing up cuttings for physical and chemical analysis. One objective of the mission was to determine the temperature of the subsurface regolith around the drill sites. Unfortunately, drilling creates a lot of heat, and experiments showed that it takes hours or even days for the soil to cool back to the original temperature. The mission’s timeline cannot afford for the rover to sit so long in one location waiting to take a measurement. Modeling may be able to help solve this problem as well. The cooling rate around the drill bit should depend on the boundary conditions, which in cylindrical coordinates centered on the drill is the ice temperature asymptotically far from the drill. If the natural subsurface temperatures are relatively constant over distances comparable to the radius that was heated by drilling, then the asymptotic temperature will equal the original temperature at the drilling location. Therefore, measuring only the cooling rate at one or several depths down the drill bit should be adequate for modeling to determine the original temperature of the subsurface. The model will need to be populated with information obtained from the drill cuttings, including the density of the soil and ice content as a function of depth, including possibly chemistry of the ice as measured by the rover’s instruments. The measured drill torque may contribute to calculating the original density of the regolith with depth. If the model has accurate constitutive equations, then with these measurements as inputs, the model can be run repeatedly using different boundary conditions until it correctly reproduces the measured cooling rates around the drill bit. This concept needs to be developed through modeling, which requires improving the model’s constitutive equations, followed by comparison to ground testing before the mission.

Mitchell and De Pater (1994) developed a one-dimensional (1D) model of heat transfer for Mercury and the Moon. It included solar insolation at the surface, geophysical heat flux from the subsurface, and constitutive equations for heat capacity and thermal conductivity of the regolith. The model used a finite-difference framework with the Crank–Nicolson algorithm. Vasavada et al. (1999) extended the model. Vasavada et al. (2012) compared the model to radiometer data of the Moon’s surface heating and cooling throughout a lunar day. Hayne et al. (2017) used it to map apparent looseness of the lunar soil globally.

Here, the model methodology is extended in three ways. This extends the progress first reported by Metzger (2018). First, the model is converted into the axisymmetric 2D Crank–Nicolson form. Second, the constitutive equations are extended based upon additional data sets for soil and mixed composition ice over varying temperatures, porosities, and gas pore pressures. Third, the heat transfer model is merged with algorithms for gas diffusion following the methodology of Scott and Ko (1968). Another model of heat and mass transfer for extraction volatiles from regolith was recently developed by Reiss (2018) using a different methodology than the one that is followed here, so comparing the two models in future work will provide a useful test of the methodologies.

The 1D thermal model described above has been reproduced and extended to the 2D axisymmetric form. The 2D heat flux equation in Cartesian coordinates iswhere $T$ = temperature; $k$ = thermal conductivity of the material; $\rho$ = density of the material; $C$ = heat capacity of the material; and $t$ = time. The equation is discretized for use in a finite-difference model. The left-hand side of the discretized equation calculates the change in $T$ from before to after one time step. The right-hand side could, therefore, be evaluated either before or after that time step. The Crank–Nicolson method is simply to average these two approaches (Crank and Nicolson 1947). This results in a linear system of equations that is stable and can be solved quickly. Using Crank–Nicolson discretization in cartesian coordinates with $\mathrm{\Delta }y=\mathrm{\Delta }x$, Eq. (1) becomes

Converting to axisymmetric form requires the extra terms in the radial derivative, so in cylindrical coordinates with $\mathrm{\Delta }r=\mathrm{\Delta }z$, it becomeswhere the radial and vertical directions are $r$ and $z$, respectively, and the discretized radial distance is $r_j=j\mathrm{\Delta }r=j\mathrm{\Delta }z$. Collecting terms with $\alpha =k\mathrm{\Delta }t/[2{(\mathrm{\Delta }z)}^2\rho C]$

Adapting the method of Summers (2012) to the axisymmetric case, two operators are defined asandso the equation becomeswhere the indices for $\alpha$ were linearized for solvability by keeping them at $n$ instead of $n+1$. The fourth-order cross-derivatives are assumed to be very small and change slowly in time relative to $\mathrm{\Delta }t$, which should be valid in realistic cases since the heat equation is diffusive and dissipative (DuChateau and Zachmann 2002)

Subtracting the left-hand side of Eq. (8) from Eq. (7) and collecting terms$T_{ij}^{*}$ is defined apart from the constants of integration by the relationshipwhich is substituted into the right-hand side of Eq. (9)

The derivatives commute$(1+{\delta }_r^2)T_{ij}^{n+1}$ and $(1-{\delta }_z^2)T_{ij}^{*}$ must therefore be equal with the correct choice of constants of integration for $T_{ij}^{*}$

Each term in this system of equationscan be represented as a tridiagonal matrix, so the tridiagonal matrix algorithm can be used to solve it efficiently.

Since this is cylindrical coordinates, the centerline $j=0$ is a special case that can be handled using the method of discretization by Scott and Ko (1968).

The model is parameterized for the thermal properties of the regolith. It also incorporates radiative heat transfer at its surface using albedo, emissivity, and insolation parameters following Mitchell and De Pater (1994) and Vasavada et al. (1999).

Parameterizing the model’s soil properties relies upon published measurements in the literature, and these measurements are cited in the following paragraphs. Those measurements show that thermal conductivity is a function of temperature, porosity (equivalently, bulk density) of the granular material, and interstitial gas pressure. The Appendix summarizes the data sets that are used in this effort.

Bulk density $\rho$ of the regolith varies on the Moon with location and depth beneath the surface, but it is not well known for asteroids. Lunar values are chosen consistent with Apollo core tubes and other Apollo measurements. Asteroid bulk densities are typically determined by fitting the results of thermal modeling to the observed thermal inertias of the asteroids. Measurements by Rosetta during a flyby of 21 Lutetia indicates the thermal inertia increases below the top few centimeters “in a manner very similar to that of Earth’s Moon” (Keihm et al. 2012). This could indicate particle sizing and/or bulk density variations over that depth, but apart from this, very little is known of possible vertical structure in asteroid regolith. The parameters in this model can be adjusted to match future spacecraft measurements to help determine asteroid regolith structure.

An important question is whether thermal conductivity also varies with particle size distribution. Chen (2008) measured thermal conductivity in terrestrial sands with different particle size distributions, varying porosity and moisture content (only the cases with zero liquid moisture are relevant to airless bodies) at ambient pressure and temperature. The $D_{50}$ median particle size of these samples varied by a factor of about 6, and samples included some with narrower (well sorted, or uniform) and broader (poorly sorted, or well graded) distributions. The results found thermal conductivity to vary with porosity but not with particle size distribution. Presley and Christensen (1997) measured thermal conductivity in soda lime borosilicate glass beads of various sizes, varying pore gas pressure from 0.5 to 100 torr at ambient temperature. Only one porosity case was measured for each grain size, with finer particles generally forming more porous packings. Since Chen’s results showed thermal conductivity is independent of particle size, the differences in thermal conductivity measured by Presley and Christensen might actually be due to the samples’ porosities, not due to their grain sizes. In contrast, Chen used realistic geomaterials, while Presley and Christensen used spherical beads. The size of contact patches between spherical particles may be correlated to grain diameter, so there may be a grain size dependence that exists in Presley and Christensen’s data that does not exist in realistic regolith. A literature review found no measurements that varied temperatures for different particle sizes while keeping constant porosities or that varied temperatures for different porosities while keeping constant particle sizes; so for now, the results by Chen are the only guidance, and they indicate thermal conductivity does not vary with particle size as an independent variable for realistic geomaterials.

Thermal conductivity for actual lunar soil has been found to follow the formwhere K = in kelvins (temperature units); and where $\chi$, $A$, and $B$ = model parameters. For example, Apollo 12 soil sample number 12001,19 was measured by Cremers and Birkebak (1971) and is shown in Fig. 1, where the dashed line is our fit using $A=0.887\mathrm{mW}/\mathrm{m}/\mathrm{K}$ and $\chi =1.56$ ($R^2=0.9929$).

Apollo 14 soil sample 14163,133 was measured by Cremers (1972) at two different bulk densities, as shown in Fig. 2: $\mathrm{1,100}\mathrm{kg}/{\mathrm{m}}^3$ (black dots with dashed curve fit, $R^2=0.9934$) and $\mathrm{1,300}\mathrm{kg}/{\mathrm{m}}^3$ (open circles with gray curve fit, $R^2=0.9926$). These densities are only 17% different and consider the difficulty of maintaining local density in an experimental apparatus, but this appears inadequate to identify a trend.

A greater variation of densities was measured by Fountain and West (1970) using crushed basalt in 10 torr vacuum at six different bulk densities, as shown in Fig. 3: ${\rho }_1=790$ (pluses), ${\rho }_2=880$ (circle), ${\rho }_3=980$ (down-triangles), ${\rho }_4=\mathrm{1,130}$ (squares), ${\rho }_5=\mathrm{1,300}$ (diamonds), and ${\rho }_6=\mathrm{1,600}\mathrm{kg}/{\mathrm{m}}^3$ (up-triangles). Note the 980 and $880\mathrm{kg}/{\mathrm{m}}^3$ samples do not follow the trend of decreasing thermal conductivity shown by the other samples probably due to experimental uncertainty.

Many functional forms were analyzed to fit all these data into one overall equation. Two forms were found to provide excellent fit, and they are discussed in the following sections. The first is a power law of the porosities, and the second is an exponential of the porosities.

For asteroids, the volatile molecules are bound in the crystalline structure of the hydrated minerals and not generally in the form of physical ice. For the Moon, the volatiles include adsorbed molecules on the surfaces of the grains as well as solid ice mixed in the regolith. The contribution of ice to thermophysical properties of regolith is determined by its chemistry and its physical state: e.g., amorphous or crystalline, snow mixed in the pore spaces, and solid ice cobbles like hail. The thermophysical properties of amorphous ice can vary by orders of magnitude depending on density and microstructure (Mastrapa et al. 2013). Amorphous ice crystallizes exothermically when there is adequate activation energy. This has been considered a mechanism for comet outbursts as heat diffuses into the interior reaching amorphous material (Sekanina 2009). In the lunar case, impact gardening could provide the activation energy, crystallizing the deposits as it matures. For now, this thermal model is based on the geological picture presented by Hurley et al. (2012) that the ice began as a homogeneous sheet and was fragmented by impact gardening, mixing grains of pure crystalline ice among grains of otherwise dry soil. The LCROSS impact did detect crystalline ice in the ejecta (Anthony Colaprete, personal communication, 2016). If the fragments are smaller than a volume element in the model, then a volumetric mixing model is adequate. These mixing models have been investigated for icy regolith by Siegler et al. (2012).

The composition of lunar ice was calculated by Tony Colaprete (personal communication, 2016) of NASA based on LCROSS impact ejecta measurements by combining measurements from the two instruments (Colaprete et al. 2010; Gladstone et al. 2010). The calculated volatile concentrations are shown in Table 1. Hydrogen gas was detected, but it should not be stable even at the temperatures of the lunar polar craters. The hydrogen gas and hydroxyl may have been products of chemistry driven by heat of the LCROSS spacecraft impact. It is beyond our present scope to back-calculate what chemicals must have been present in the ice prior to the impact. For now, this remains the best estimate of the composition of lunar ice.

The saturation curves of these volatiles (NIST 2017) shown in Fig. 9 illustrate that temperatures adequate to release water from the regolith will also release many other volatiles. For now, only water sublimation has been modeled. Modeling here may treat the sublimation of each species separately based on partial pressures and treat the diffusion of gas through the pore spaces based on overall pressure and molecular collision rates in the mixed gas. This assumption needs to be checked with measurements of actual lunar ice. Kouchi et al. (2016) found that ice mixtures of CO to ${\mathrm{H}}_2\mathrm{O}$ in ratios $50:1$ and $10:1$ subjected to conditions for the sublimation of the CO but not the ${\mathrm{H}}_2\mathrm{O}$ left the water ice in a porous amorphous state with a density similar to high-density amorphous ice. They also found that at 140 K this matrix-sublimated high-density amorphous ice transitioned to cubic ice. Doubtless, the porosity resulting from matrix sublimation decreases thermal conductivity of the remaining matrix and transitioning back to crystalline ice increases it, so these effects may occur when subliming mixed composition lunar ice. The nonwater species of lunar ice constitute less than 50% by weight of the combined ices, so the induced porosity should be much less than reported by Kouchi et al. (2016). For now, the effect is ignored. Future work may add an ad hoc parameter to treat it simplistically, but it would be little more than a guess. To inform a better model, experimental measurement is needed for the thermal conductivity of matrix sublimed, mixed composition ice.

For thermal conductivity, the contribution of pure crystalline water ice is calculated using the data points of Ehrlich et al. (2015), reproduced in Fig. 10 with the added curve fitin $\mathrm{W}/\mathrm{m}/\mathrm{K}$, and where temperature $T$ is in kelvins. Available thermal conductivity data for the other volatiles are limited. Sumarokov et al. (2009) measured CO ice in the range 1–20 K. At 20 K, it is about $0.5\mathrm{W}/\mathrm{m}/\mathrm{K}$, more than an order of magnitude less than the extrapolation of water to that temperature per Fig. 10. Koloskova et al. (1974), cited in Sumarokov et al. (2003), measured ${\mathrm{CO}}_2$ and found it about $1\mathrm{W}/\mathrm{m}/\mathrm{K}$ at 100 K, about 1/6 the value of water. Manzhelii et al. (1972) report solid ammonia about $1.8\mathrm{W}/\mathrm{m}/\mathrm{K}$ at 100 K, about 1/3 the value of water. Lorenz and Shandera (2001) found ammonia-rich (approximately 10%–30%) water ice has a thermal conductivity of about 1/2 to 1/3 that of pure water ice.

In the 1D model of sandwiched materials, the net thermal conductivity iswhere $d_i$ and $k_i$ = thickness and thermal conductivity of each layer. To first-order approximation, which is the limit of accuracy considering the other unknowns, a mixture of dry regolith with ice grains would scale asin $\mathrm{W}/\mathrm{m}/\mathrm{K}$, where $\sim 23\%$ volume fraction of ice is derived from $\sim 8.9\%$ by weight of ice, a rough estimate based on Table 1 approximating the mixed chemistry as if it were all water. This indicates $k_{\mathrm{bulk}}\sim 0.0013$, only about 30% higher than dry regolith. If $k_{\mathrm{ice}\mathrm{grains}}=2\mathrm{W}/\mathrm{m}/\mathrm{K}$ to reflect the mixed chemistry instead of $6\mathrm{W}/\mathrm{m}/\mathrm{K}$ for pure water ice, then $k_{\mathrm{bulk}}\sim 0.0013$, not measurably changed. The mixed chemistry of ice cobbles can safely be ignored for modeling thermal conductivity.

As volatiles are released, the increasing pore pressure increases the thermal conductivity by orders of magnitude, as shown in Fig. 8, where FW is in hard vacuum, PC is in a range of pore pressures, and Chen is at Earth ambient pressure. The order of magnitude of $k$ compares reasonably for all data sets when pore pressure is accounted for, although the shape of PC disagrees with the other two, as discussed previously. There are inadequate data regarding how to reconcile this, but the following observations lead to a hypothesis. First, as shown in Fig. 8, the PC data with higher porosity are better distributed between the end points formed by the Chen and FW curves than they are at the lower porosities. Second, as shown in Fig. 11, the data at high porosity are continuous with the upper pressure end-point, while the low porosity data are discontinuous. Fig. 11 plots the ratiowhere $k_{\mathrm{PC}}(\nu ,P)$ = PC data; and $k(\nu ,300\mathrm{K})$ = Eq. (26) based on FW data evaluated at $T=300\mathrm{K}$ (the temperature chosen to match the temperature of the PC data). $k_{\mathrm{PC}}$ is evaluated at $P=0.5\mathrm{Pa}$ in the denominator, which is apparently below the floor where pore gas does not contribute significantly to thermal conduction in the soil, as discussed in the following section. Appended on the right side of each set of points is one data pointwhere $k_{\mathrm{Chen}}(\nu )$ = Eq. (22) with the fitting parameters found by Chen and zero moisture content. The continuity for the high porosity cases suggests a hypothesis that the high porosity cases are correct, while the low porosity cases (coarse particles with less cohesion) suffered experimental disturbance. Third, it is understandable why the more porous cases of PC would be more accurate than the less porous cases because they are more stabilized by higher cohesion, as discussed.

Consistent with the assumption that the most porous data of PC are the least disturbed in the laboratory measurements, trendlines were projected on log-log axes in Fig. 12 such that they intersect at the same point where the FW and Chen trendlines are projected to intersect. These projections pass through the most porous case of PC data. Dots on the top and bottom trendlines are calculations using the Chen and FW fitted functions at the porosities of the PC data.

The assumption is that these trend lines are what would have been measured in PC had there been no mechanical disturbance of the samples. This assumption is necessary to reconcile the existing data sets and create a constitutive equation. The family of trend lines is viewed in Fig. 13 through two different projections: into the ($\nu ,k$) plane and into the ($P,k$) plane, where $P$ is pore pressure in the soil. In the ($\nu ,k$) plane, the top line is for 760 torr pore pressure (Chen), the bottom is for ${10}^{-8}\mathrm{torr}$ (FW), and the intermediate are for 100 torr (upper) to 0.5 torr (lower) (PC). These trendlines were chosen to intercept at the same point off the right side of the figure where FW and Chen also intercept.

In the ($P,k$) plane, the far right vertical column of points is from Chen, the left vertical column of points is from FW, and the second column of points from the left is not from any dataset but is the point where the PC data project to an intersection with their corresponding floor. Each floor represents conduction and radiation through the grains without a significant contribution from pore gas. Note that the existence of this floor implies that the gas contribution is additive to the other contributions (not multiplicative).

The general fitting function for this family of curves iswhere $\hat{P}=\mathrm{Max}(P,P_0)$; and $P_0$ = pressure below which is the floor of minimum conductivity. At constant pressure, this reduces to the form of a single exponentialwhereas Eq. (26) at constant temperature reduces to the form of a double exponential

The second term is the coefficient for the radiation term in $T^3$. Radiation ought to be independent of gas pressure to good approximation in the rarefied conditions considered here, so the additive form of Eq. (34) agrees with expectations. Eq. (32) was developed from data measured at $T=300\mathrm{K}$. Therefore, the $T^3$ in Eq. (26) can be added to these curve-fits after subtracting the assumed ${(300\mathrm{K})}^3$ contribution. With some manipulation, this yields the full modelin $\mathrm{W}/\mathrm{m}/\mathrm{K}$, where curve fitting with $T$ in kelvins and $P$ in pascals provided the following constants:

The decimal places are not all significant, but they are the exact values coded into the model. Quantifying significant digits of model parameters is left to future work when better empirical datasets are available and when the knowledge gaps in the physics have been reduced.

For specific heat, this model uses a mass-weighted mixing model of ice and regolith. The contribution of the dry regolith is informed by the measurements previously made for Apollo soil samples and analog materials. Fig. 14 shows a representative comparison. Apollo samples 10084 and 10057 are from Winter and Saari (1969). The empirical fitting functions by Winter and Saari (1969)and by Hemingway, Robie, and Wilson (1973)both in $\mathrm{J}/\mathrm{kg}/\mathrm{K}$, are very close to one another and only slightly higher than the experimental data above 250 K. New fifth-order and fourth-order polynomial fits were tried. The fifth order diverges from the probable trend just outside the range of experimental measurements, so it is rejected. The fourth order is marginally better than the one by Hemingway, Robie, and Wilson (HRW) and seems to preserve the trends in extrapolation. The fit by HRW is actually based on a larger set of measurements, including Apollo soil samples 14163, 15301, 60601, and 10084 to represent the average of lunar soil, so HRW is selected. It fits the data with less than a 10% error.

The heat capacity of dry regolith should vary proportionally to the soil’s solid fraction ($1-\nu$), which varies by about $\pm 16\%$ from the median value in the lunar case. For the asteroid case, the bulk density may vary widely due to large changes in particle size and ultra-low gravity. Measurement of asteroid regolith density in situ, including any stratigraphic variation in the subsurface, is required to guide more accurate models.

The specific heat of pure water ice was measured by Giauque and Stout (1936) from about 15 to 270 K and isin $\mathrm{J}/\mathrm{kg}/\mathrm{K}$, where T is in kelvins.

The specific heats of the major volatiles in lunar ice are shown in Fig. 15: (in order of prevalence) water by Giauque and Stout (1936), hydrogen sulfide by Giauque and Blue (1936), sulfur dioxide by Giauque and Stephenson (1938), ammonia by Overstreet and Giauque (1937), carbon dioxide by Giauque and Egan (1937), ethylene by Egan and Kemp (1937), methanol by Carlson and Westrum (1971), methane by Colwell et al. (1963), and carbon monoxide measured by Clayton and Giauque (1932). Weighting these according to Table 1, the composite heat capacity is shown in Fig. 16. This neglects the hydrogen and hydroxyl that were also measured in the lunar ice ejecta, which are assumed to have come from the decomposition of unidentified components.

The composite heat capacity for ice is calculated by a mass-weighted sum of the individual components. This assumes a linear mixing model which may not be correct depending on the actual crystalline or amorphous form of the ice, but until measurements are taken on the Moon, this is the best assumption that can be made. Above the sublimation temperature of each component, the weighting is renormalized for the reduced mass. The heat capacity for the composite ice and for pure water are compared in Fig. 16. The integrated heat capacity of pure water is found to be always within 14% of the integrated heat capacity for composite ice. In the present accuracy of the approximation, pure water’s heat capacity can be used as an adequate representation of the composite ice. The combined specific heat of regolith with 8.9% by weight ice (now approximating it is all water) is shown in Fig. 17. The specific heat of water ice is roughly three times higher than the specific heat of dry lunar soil, but since it constitutes only 8.9% by weight of the regolith, it raises total heat capacity of the mixture by only about 29% at 40 K and about 16% at 200 K.

The sublimation of water ice on the Moon at temperatures below the triple point is treated by Andreas (2007) by relating the saturation vapor pressure of water ice, $e_{\mathrm{sat},\mathrm{i}}(T)$, which is a function of temperature $T$, to the sublimation rate $S_0$in $\mathrm{kg}/{\mathrm{m}}^2/\mathrm{s}$, where $M_W$ = molecular weight of water; and $R$ = universal gas constant. Kossiacki and Leliwa-Kopystynski (2015) modified this by subtracting the partial pressure of the vapor $P$ from the saturation pressure

When partial pressure reaches saturation pressure, then sublimation should cease. Here, the model uses the vapor pressure relationship provided by Murphy and Koop (2005)in pascals.

The free surface area of the ice where sublimation takes place depends on the physical state of the ice, whether it exists as large cobbles of ice surrounded by regolith fines, fine particles of ice comparable to the size of regolith particles intermixed with the mineral grains, a rind of ice coating the mineral grains, amorphous material residing in the pore spaces between grains, or another form. How it is modeled depends on which physical state is assumed. One simple way to model the exposed surface area of the ice $\mathrm{\Delta }s$ (in a cell toroidal of radius $r$ about the model’s axis because this is an axisymmetric model) is $\mathrm{\Delta }s=2\pi r\mathrm{\Delta }z\sigma$, where $\sigma$ = parameter based on the expected physical state of the ice informed by the geological model. The net sublimed mass during the $n$th timestep in cell location ($i,j$) is therefore

The mass of water in the regolith’s pore spaces in the toroidal cell could be modeled aswhere ${\phi }_{ij}$ = fraction of the pore space filled by ice; and ${\rho }_I$ = density of the ice. A relationship is needed between ${\sigma }_{ij}$ and ${\phi }_{ij}$ to represent the physical state of the ice. For now, the model uses the simplification ${\sigma }_{ij}={\phi }_{ij}$. Heat capacity is a simple scaling between ice and regolith heat capacities by the amount of each mass within a cell. The temperature in a cell may continue to rise even as ice sublimes until it reaches the triple point because sublimation is a slow process at these temperatures and the system remains in nonequilibrium. In practical cases that have been modeled, the temperature never rose as high as the triple point. As sublimation occurs, the heat of fusion is subtracted from the internal energy of the cell, and the temperature is lowered accordingly before time-stepping the model to calculate the conduction of heat again. The gas and the solid components in each cell are assumed to have the same temperature.

Hydrated minerals in asteroid regolith release their volatiles as a function of temperature. For testing an asteroid mining prototype, it was economically beneficial to use lower temperature materials for early tests, so a lower temperature simulant was developed using primarily epsomite because it releases most of its water of hydration below 300°C. Curves were obtained through thermo-gravimetric analysis (TGA) to determine the mass of released volatiles at each temperature increment. Examples of these curves are shown in Fig. 18 for asteroid simulant UCF-CI-1 measured by Metzger et al. (2019), the Orgueil meteorite by King et al. (2015), and epsomite by Ruiz-Agudo et al. (2007).

To model this, as a region in the model reaches a new high temperature, the volatiles up to that temperature per the TGA curve are released as vapor, and the model remembers that no more volatiles will be released from that location until an even higher temperature is achieved. The appropriate amount of gas per the TGA curve is added to the gas already in the pore space at that location. Energy spent liberating the volatiles in each step is removed from the regolith appropriately in each time step. Thermal and gas diffusion are then iterated. The equation to fit epsomite [Fig. 18(b)] iswhere $w$ = weight percent (of a cell’s material) that has sublimed; and $T$ = in kelvins.

The model does not modify the thermal conductivity of the solid fraction of the soil as mass is converted to vapor, although it should reduce that term because (especially with epsomite) a large fraction of the solid mass is lost, reducing the solid conduction contact network. That effect is offset by the increase in conductivity due to rising pore pressure as shown in Fig. 13, but there are no empirical data at present to guide this improvement. Those experimental measurements and modeling are left to future work.

The model incorporates diffusion using the finite-difference equations of Scott and Ko (1968). As vapor is evolved as described, the pressure differences drive it into neighboring cells. To couple the fast gas diffusion equations and the slow thermal diffusion equations while maintaining the stability of the model, it was necessary to implement different time steps for each set of equations. Adaptive timesteps were implemented for the fast diffusion process, dividing each heat flow timestep into the minimum number of smaller timesteps necessary for a stable solution of the gas diffusion equations. As pressure gradients increase, the gas diffusion timesteps become smaller. The resulting model is fast, allowing simulation of an hour-long physical test in just 5 or 10 min on a standard laptop computer.

Only limited testing of the model has been performed. The following four cases demonstrate aspects of the thermal algorithms and the overall code structure with increasing complexity. The first case is 1D simulations of the lunar regolith heating and cooling in sunlight at various latitudes in comparison with Lunar Reconnaissance Orbiter (LRO) Diviner data per Fig. 9(a) of Vasavada et al. (2012). The lunar surface albedo as a function of angle and other parameters choices by Vasavada et al. (2012) were intertwined with choices of thermal conductivity to make their model match lunar data sets. Vasavada et al. (2012) used $k=0.6\mathrm{mW}/\mathrm{m}/\mathrm{K}$ for the most porous soil at the lunar surface (${\nu }_0=0.58$, bulk density ${\rho }_0=\mathrm{1,300}\mathrm{kg}/{\mathrm{m}}^3$), asymptotically approaching $k=7\mathrm{mW}/\mathrm{m}/\mathrm{K}$ for the least porous soil at depth (${\nu }_{\infty }=0.42$, bulk density ${\rho }_{\infty }=\mathrm{1,800}\mathrm{kg}/{\mathrm{m}}^3$), with the porosity exponentially decaying as a function of depth, $z$where $H$ = single remaining model parameter. Values of $H$ can be iterated until the model makes predictions that match observations of lunar surface temperature rising and falling as the Moon rotates in the sunlight. The value of $H$ is thus a proxy to characterize how rapidly the soil compactifies with depth at each location on the Moon in some averaged sense.

Following the choices of Vasavada et al. (2012), the thermal conductivity form in Eq. (15) becomes

Ice content is set to zero. The specific heat is represented by Eq. (38). Simulations were performed for the Moon rotating in the sunlight over 37 lunations (months) to achieve a steady state. The final lunation is shown in Fig. 19 for three cases: $H=0.5\mathrm{cm}$ (short dots, top curve, most compacted soil so highest thermal inertia), 3.5 cm (best fit, solid curve), and 30 cm (long dashes, bottom curve, loosest soil so least thermal inertia). The model was successful in predicting lunar temperatures, indicating the model is structured correctly. Future work will use the improved parameterization of Eq. (25), which will make it necessary to determine how albedo and the other model parameters must be changed from the values of Vasavada et al. (2012) to match lunar measurements. This should produce improved characterization of $H$.

The second case is for equatorial conditions on asteroid 101955 rotating in sunlight. Fine-tuning of the model has not been performed for the asteroid because adequate data sets from asteroids do not exist, but better data are expected soon from spacecraft missions currently in progress. Parameterization is therefore speculative. Many cases were modeled, and they produced similar results with differences that can be tested when the mission data become available. This particular case shown in Fig. 20 used a three-layer regolith model assuming the surface and deepest layers of the asteroid have identical properties, while a layer with different properties exists between 0.5 and 6.0 cm depth. Theory says such an intermediate layer might form on asteroids by thermal cracking as the asteroid rotates in the sunlight, while the uppermost layer loses the fines in the low gravity. The surface and deepest layers were assumed to be very porous with bulk density $\rho =\mathrm{1,300}\mathrm{kg}/{\mathrm{m}}^3$, while the intermediate layer was assumed to have $\rho =\mathrm{2,340}\mathrm{kg}/{\mathrm{m}}^3$. The specific heat was assumed the same as lunar soil in all layers per Eq. (38). The thermal conductivities of all three layers were assumed to follow Eq. (15). For the surface and deepest layers, parameter $A$ was estimated by taking the value of Vasavada et al. (2012) for the most porous lunar soil, $k=0.6\mathrm{mW}/\mathrm{m}/\mathrm{K}$, and then multiplying by the particle size factor suggested by Presley and Christensen (1997), ${(D_{\mathrm{asteroid}}/D_{\mathrm{lunar}})}^{0.5}$, where $D_{\mathrm{asteroid}}$ = average particle size of the asteroid regolith; and $D_{\mathrm{lunar}}$ = average particle size of lunar soil. This could be interpreted as the expected larger contact patches because asteroid regolith is dominated by large gravel particles. $D_{\mathrm{asteroid}}=1.5\mathrm{cm}$ and $D_{\mathrm{lunar}}=60\mathrm{\mu m}$ result in $A=9.49\mathrm{mW}/\mathrm{m}/\mathrm{K}$. $\chi =2.7$ is kept, matching Vasavada et al. (2012). The thermal conductivity of the intermediate layer was assumed to have $A=3.736\mathrm{mW}/\mathrm{m}/\mathrm{K}$, corresponding to an intermediate value between the most and least compacted lunar soil, and $\chi =0.434$ for reduced radiative heat transfer due to reduced porosity. The simulation replicated the solar insolation conditions for Bennu and its approximately 4.3 h rotation rate for 500 rotations to achieve steady state. The resulting range of temperatures shown in Fig. 20 correctly matches the range observed on Bennu, as it rotates in the sun per Lauretta et al. (2015).

The third case adds complexity by using the 2D axisymmetric version of the Crank–Nicolson formulation while retaining the lunar soil property equations of the 1D model [following Vasavada et al. (2012)]. It is a simulation for a drilling test in which a warm drill bit is embedded in soil that carries away its heat. The soil is in a tall, narrow, cylindrical container 14 cm in radius and 120 cm tall. The experiment is inside a warm vacuum chamber. This is the geometry of simulated tests done with a Honeybee Robotics drill at a NASA Glenn Research Center vacuum chamber where a liquid nitrogen bath kept the soil container at constant temperature (77 K) and removed heat from the soil conductively.

In these simulations, four different boundary temperatures (133, 153, 173, and 193 K) instead of the liquid nitrogen bath temperature were successively used to test how the Resource Prospector Mission drill bit could measure cooling rate while embedded in soil. The initial soil temperature before drilling was set to the boundary temperature. The soil model had no ice and was set to ${\rho }_0=\mathrm{1,300}\mathrm{k}/{\mathrm{m}}^3$ (${\nu }_0=0.58$), ${\rho }_{\infty }=\mathrm{1,950}\mathrm{k}/{\mathrm{m}}^3$ (${\nu }_{\infty }=0.37$), $H=5\mathrm{cm}$, porosity following Eq. (46), heat capacity per Eq. (38), thermal conductivity per Eq. (15) parameterized by

Albedo and emission at the surface follow Vasavada et al. (2012). The vacuum chamber walls are set to 193 K for radiative heat transfer. The drill bit was initially held at 213 K for 5,000 time steps (2 s each) while the soil came to equilibrium, shown in Fig. 21. Then, it was allowed to cool to determine the cooling rate of the bit. The bit and drilling mechanism attached to its upper end were modeled for realistic thermal inertia. The model showed that the cooling process is so slow in lunar soil that it takes hours for soil around a warm drill bit to return to ambient temperature, matching the experimental observations. It is impractical for a lunar rover to pause its mission until the soil cools to take a subsurface temperature measurement. Fig. 22 shows that the cooling rate of the bit depends on the boundary temperature condition, which represents the original subsurface soil temperature before drilling. Thus, the lunar drill can quickly measure the cooling rate and continue its mission, relying on modeling to calculate the original subsurface temperature. The actual lunar case would be more complex than what was simulated here because boundary temperatures (temperature asymptotically far from the drill) should vary with depth. Additional uncertainty exists from the compaction of the soil and ice content with depth. These additional unknowns will be informed in part by drill torque as it is inserted into the subsurface and by analysis of the cuttings for ice content as the drill brings the cuttings up to instruments at the surface. All these measurements would need to be analyzed to inform model parameterization. A future study could analyze the accuracy of this overall method and its sensitivity on the several parameters.

The fourth case integrated the full set of new constitutive equations given by Eqs. (36), (38), and (45) into the 2D axisymmetric Crank–Nicolson model. This was used to simulate the extraction of water from asteroid regolith by the WINE spacecraft coring device (Zacny et al. 2016; Metzger et al. 2016). The corer is a hollow tube with flutes on the exterior that enable it to drill into the subsurface, filling the hollow center of the tube with regolith. The corer walls are heated on the inside, while its layered insulating structure minimizes heat transfer to its exterior. Regolith increases in thermal conductivity as it warms, so the process becomes increasingly efficient. It releases volatiles according to the TGA curve, further increasing thermal conductivity.

The simulations were performed both for terrestrial test conditions with a 1 bar background pressure in the regolith and for space applications with the pores initially in vacuum. Fig. 23 is a case inside of a vacuum chamber with zero ambient pressure and wall temperatures at 293 K. The soil begins at 272 K temperature, and the soil container’s boundaries are kept at 272 K throughout the simulation. In each timestep, 200 W thermal energy is delivered to the inside walls of the corer uniformly along its heated surface and then diffused from the tube through the soil.

Videos were created from the simulation data showing the resulting temperature and pressure fields in the regolith. Fig. 23 shows a series of snapshots of the temperature field in the cross section through the corer. The simulation demonstrated that corer design successfully keeps most of the thermal energy inside the interior, although some energy leaks to the exterior. The videos show that the pressure builds up almost immediately and then decays and the pressure field becomes more uniform. This decay is because the simulant’s water of hydration becomes depleted. Fig. 24 shows a series of snapshots of the vapor pressure field. The semicircular pressure gradient at the top inside the corer is where the vapor diffuses to the collection tube located on the centerline. In the corresponding experiments, the tube leads to the cold trap where volatiles are frozen, keeping the tube at near vacuum conditions, but in the simulation, the vapor that reaches the tube’s entrance is simply accounted for and then removed from the simulation to maintain the tube entrance at vacuum conditions. In the initial simulations, a significant fraction of the vapor can be seen exiting the bottom of the corer rather than diffusing into the collection tube, reducing the system’s mining efficiency. In this particular case, the soil’s initial temperature started near the triple point, so as it was warmed, the vapor exiting the bottom of the corer did not freeze elsewhere in the soil but diffused to the soil’s upper surface where it escaped into the surrounding vacuum.

To study how to capture a larger fraction of the vapor, additional simulations were performed where a gap was left between the top of the soil and the inside top of the coring tube. This can be achieved experimentally by not driving to the corer all the way into the soil. This is the case shown in Fig. 25. The gap is so small it is not visible, but it is simulated by appropriate choice of model parameters so the gas can diffuse out from the soil into vacuum all along its top surface inside the corer. This produced a flat pressure gradient across the entire top of the soil inside the corer instead of the semicircular pressure gradient in Fig. 24. Comparing with Fig. 24, the vapor pressure outside the coring tube was reduced because vapor was transported upward through the corer more efficiently. This increased water capture by 520%. This illustrates how modeling can be used to drive the design of mining devices for improved performance in an extraterrestrial environment.

Thermal volatile extraction modeling has been successfully developed for the 1D and axisymmetric 2D cases for asteroid and lunar regolith. The modeling includes parameterization for regolith thermal conductivity and heat capacity based on measurements of lunar soil samples, simulants, and terrestrial soil and ices. This is apparently the first time a soil constitutive model has successfully reconciled datasets for temperature, porosity, and gas pore pressure variables into a single equation. The model has been only partially tested. It produced excellent agreement with LRO Diviner data of the Moon and estimates of asteroid Bennu heating and cooling as they rotate in the sun. The 2D axisymmetric features have been demonstrated by simulating the Resource Prospector drill cooling after insertion into lunar soil. The model can also simulate the effects of ices upon the thermal conductivity and heat capacity of the lunar regolith. (For the asteroid case, ice is not expected as the volatiles are in the form of hydrated minerals.) The model is based on the assumption that lunar ice is crystalline rather than amorphous, which is supported by some data although the presence of amorphous phases cannot be ruled out. More work is needed to adapt the model to include the effects of amorphous ice. The model also successfully integrated equations for volatile release and gas diffusion along with the thermal diffusion equations, employing multiple adaptive time steps to handle the different characteristic times of each part of the physics. The fully integrated model has been demonstrated for the case of a corer heating and extracting volatiles from asteroid regolith.

The following symbols are used in this paper:

$A$, $B$

model fitting coefficients defined in the text;

$a$, $b$

model fitting exponents defined in the text;

$\hat{a},\hat{b},\hat{c}$

model fitting parameters defined in the text;

$C$

soil heat capacity;

$c_1,\dots ,c_7$,

model fitting exponents defined in the text;

$d$

soil layer thickness;

$d_2,\dots ,d_4$

model fitting exponents defined in the text;

Exp

exponential function;

$e_{\mathrm{sat},\mathrm{i}}$

saturation vapor pressure of water ice;

$i$

index of model cell location in vertical direction;

$j$

index of model cell location in radial direction;

$k$

thermal diffusivity constant;

$k_1,\dots ,k_9$

model fitting exponents defined in the text;

$k_{\mathrm{Chen}}$

thermal diffusivity values measured by Chen (2008);

$k_{\mathrm{PC}}$

thermal diffusivity values measured by Presley and Christensen (1997);

Ln

natural logarithm function;

$M_W$

molecular weight;

$m_{\mathrm{S}}$

mass of ice sublimed;

$n$

index of model timesteps;

$O()$

order of magnitude;

$P$

pore gas pressure in the soil;

$R$

universal gas constant;

$r$

location in radial direction;

$S_{\mathrm{r}}$

moisture saturation of soil;

$S_0$

sublimation rate of ice;

$T$

temperature;

$t$

time;

$V$

volume;

$w$

weight percent sublimed;

$z$

depth into the soil column;

$\alpha$

model thermal parameter defined in the text;

$\mathrm{\Delta }$

model step difference in variable $r$, $z$, or $t$;

${\delta }_z^2,{\delta }_r^2$

calculus operators defined in the text;

$\epsilon$

model fitting exponent defined in the text;

$\nu$

soil porosity;

$\pi$

pi (the number);

$\rho$

bulk density of soil;

${\rho }_I$

density of water ice;

$\sigma$

parameter to define the surface area of ice in a model cell;

$\phi$

fraction of soil’s mass that is ice; and

$\chi$

radiative transfer coefficient.

Some or all data, models, or code generated or used during the study are available from the corresponding author by request: Mathematica notebook containing data from Figs. 1–8, 10–14, 16–20, and 23–25; Excel spreadsheet containing data from Figs. 9 and 15.

This work was directly supported by NASA SBIR contract no. NNX15CK13P, “The World is Not Enough (WINE): Harvesting Local Resources for Eternal Exploration of Space.” This work was directly supported by NASA’s Solar System Exploration Research Virtual Institute cooperative agreement award NNA14AB05A.