Mars: Quantitative Evaluation of Crocus Melting behind Boulders

By definition, the molecular flux (index m) and convective flux (index c) are related by

$$\begin{eqnarray}&&{Q}_{c}={Q}_{m}\mathrm{Nu}\qquad {\rm{heat}}\,{\rm{flux}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{E}_{c}={E}_{m}\mathrm{Sh}\qquad {\rm{mass}}\,{\rm{flux}},\end{eqnarray} \tag{ 2 }$$

which defines the Nusselt number, Nu, and the Sherwood number, Sh. The molecular heat flux, Q_m, is proportional to the temperature difference, ΔT, and the convective heat flux, Q_c, can be written in the same form, only the coefficient of proportionality is no longer constant:

$$\begin{eqnarray}&&{Q}_{m}=\displaystyle \frac{{k}_{m}}{L}{\rm{\Delta }}T=\displaystyle \frac{{\kappa }_{m}}{L}\rho {c}_{p}{\rm{\Delta }}T,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{Q}_{c}=\displaystyle \frac{{k}_{c}}{L}{\rm{\Delta }}T=\displaystyle \frac{{\kappa }_{c}}{L}\rho {c}_{p}{\rm{\Delta }}T.\end{eqnarray} \tag{ 4 }$$

Here, L is a characteristic length scale, practically the dimension of the plate the laboratory experiment is conducted with. In some situations, L eventually cancels from the final expressions. For the mass fluxes,

$$\begin{eqnarray}&&{E}_{m}=\displaystyle \frac{{D}_{m}}{L}{\rho }_{w},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{E}_{c}=\displaystyle \frac{{D}_{c}}{L}{\rho }_{w},\end{eqnarray} \tag{ 6 }$$

where D_c is no longer constant. Here ρ_w is the density of water vapor, and the subscript w refers to H₂O.

The governing equations for heat transfer without mass exchange and mass exchange without heat flow are, respectively,

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}T={\kappa }_{m}{{\rm{\nabla }}}^{2}T,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial p}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}p={D}_{m}{{\rm{\nabla }}}^{2}p,\end{eqnarray} \tag{ 8 }$$

where p is the partial pressure. The mechanical diffusivity, D_m, is that of one substance in another, D₁₂ = D₂₁ (binary diffusion), not that of a substance within itself, D₁₁ (self diffusion). For free convection, the velocity ${\boldsymbol{u}}$ is in turn given by linear relations with T or p, respectively.

As the equations governing heat transport and mass transport are mathematically equivalent, both have the same solutions, even for turbulent flow. Both transport coefficients are given by the same function Φ (Jakob 1949),

$$\begin{eqnarray}&&{\kappa }_{c}={\kappa }_{m}{\rm{\Phi }}(\mathrm{Gr},\Pr ),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{D}_{c}={D}_{m}{\rm{\Phi }}(\mathrm{Gr},\mathrm{Sc}).\end{eqnarray} \tag{ 10 }$$

The dimensionless numbers Gr, Pr, and Sc are similarity variables and will be defined below. In other words, $\mathrm{Nu}={\rm{\Phi }}(\mathrm{Gr},\Pr )$ and $\mathrm{Sh}={\rm{\Phi }}(\mathrm{Gr},\mathrm{Sc})$. Due to this similarity relation, measurements of ${\kappa }_{c}/{\kappa }_{m}$ for heat transfer can be used to determine D_c/D_m for mass transfer.

The function Φ has to be determined from laboratory measurements. For a large horizontal plate, Holman (1990) quotes the experiments by Fujii & Imura (1972) in water:

$$\begin{eqnarray}&&\mathrm{Nu}=0.13\,{\mathrm{Ra}}^{1/3}=0.13{\left(\mathrm{Gr}\times \Pr \right)}^{1/3},\end{eqnarray} \tag{ 11 }$$

$\mathrm{Ra}=\mathrm{Gr}\,\times \,\Pr$ is the Rayleigh number. Incropera et al. (2007) list $\mathrm{Nu}=0.15\,{\mathrm{Ra}}^{1/3}$ for high Rayleigh number. Jakob (1949), which is the source of Ingersoll (1970), lists

$$\begin{eqnarray}&&\mathrm{Nu}=0.068\,{\mathrm{Gr}}^{1/3}\end{eqnarray} \tag{ 12 }$$

based on laboratory measurements of heat transfer in air (Mull & Reiher 1930).

In the turbulent regime (at high Grashof number), the universal function Φ is in all these cases of the form

$$\begin{eqnarray}&&{\rm{\Phi }}(\mathrm{Gr},N)=C\,{\mathrm{Gr}}^{1/3}{N}^{1/3},\end{eqnarray} \tag{ 13 }$$

where N can be either $\Pr$ or $\mathrm{Sc}$. Niemela et al. (2000) measured

$$\begin{eqnarray}&&\mathrm{Nu}=0.124\,{\mathrm{Ra}}^{0.309\pm 0.004}\end{eqnarray} \tag{ 14 }$$

in liquid helium and a specific container geometry. The small deviation from 1/3, if it applies to all geometries, would cause the length scale L to appear in the final result. From Equations (2), (10), and (13), we obtain the desired relation

$$\begin{eqnarray}&&\mathrm{Sh}=C{\left(\mathrm{Gr}\times \mathrm{Sc}\right)}^{1/3}.\end{eqnarray} \tag{ 15 }$$

The Grashof number is defined by

$$\begin{eqnarray}&&\mathrm{Gr}={L}^{3}\displaystyle \frac{g}{{\nu }^{2}}\displaystyle \frac{{\rm{\Delta }}\rho }{\rho },\end{eqnarray} \tag{ 16 }$$

where g = 3.71 m s⁻² is the specific surface gravity on Mars. The relative density difference, ${\rm{\Delta }}\rho /\rho$, for temperature-driven buoyancy and humidity-driven buoyancy, respectively, is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\rho }{\rho }=\beta {\rm{\Delta }}T=\displaystyle \frac{{\rm{\Delta }}T}{T},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\rho }{\rho }=\displaystyle \frac{{p}_{w}({M}_{1}-{M}_{w})}{{p}_{0}{M}_{1}-{p}_{w}({M}_{1}-{M}_{w})},\end{eqnarray} \tag{ 18 }$$

where β is the coefficient for thermal expansion, which is 1/T for an ideal gas. In Equation (17), Δρ is the density difference between the warm and the cold air, whereas in (18) it is the density difference between the humid and the dry air.

The formula for the mass flux (6) then becomes

$$\begin{eqnarray}&&{E}_{c}={{CD}}_{m}{\rho }_{w}{\left(\displaystyle \frac{g}{{\nu }^{2}}\displaystyle \frac{{\rm{\Delta }}\rho }{\rho }\right)}^{1/3}{\left(\displaystyle \frac{\nu }{{D}_{m}}\right)}^{1/3},\end{eqnarray} \tag{ 19 }$$

where the length scale L has now canceled, and $\mathrm{Sc}=\nu /{D}_{m}$. The quantity

$$\begin{eqnarray}&&\gamma ={\left(\displaystyle \frac{{\nu }^{2}}{g}\displaystyle \frac{\rho }{{\rm{\Delta }}\rho }\right)}^{1/3}\end{eqnarray} \tag{ 20 }$$

can be called "Grashof length scale" because it has units of length. One can also define $C^{\prime} =C{\mathrm{Sc}}^{1/3}$, although unlike C, $C^{\prime}$ is not a universal constant but depends on the Schmidt number for the specific gas composition. With these abbreviations,

$$\begin{eqnarray}&&{E}_{c}={D}_{m}\displaystyle \frac{C}{\gamma }{\mathrm{Sc}}^{1/3}{\rho }_{w}={D}_{m}\displaystyle \frac{C^{\prime} }{\gamma }{\rho }_{w}=\displaystyle \frac{L}{\gamma }{\mathrm{Sc}}^{1/3}{E}_{m}.\end{eqnarray} \tag{ 21 }$$

Mills (2001), used in Dundas & Byrne (2010), directly provides

$$\begin{eqnarray}&&\mathrm{Sh}=0.14\,{\left(\mathrm{Gr}\times \mathrm{Sc}\right)}^{1/3}.\end{eqnarray} \tag{ 22 }$$

Hecht (2002) uses

$$\begin{eqnarray}&&\mathrm{Sh}=0.15\,{\left(\mathrm{Gr}\times \mathrm{Sc}\right)}^{1/3}.\end{eqnarray} \tag{ 23 }$$

The expression from Ingersoll (1970) is

$$\begin{eqnarray}&&{E}_{c}=0.17{\rho }_{w}{D}_{m}{\left(\displaystyle \frac{g}{{\nu }^{2}}\displaystyle \frac{{\rm{\Delta }}\rho }{\rho }\right)}^{1/3}.\end{eqnarray} \tag{ 24 }$$

Compared with Equation (19), this amounts to

$$\begin{eqnarray}&&C\,{\mathrm{Sc}}_{{{\rm{H}}}_{2}{\rm{O}}\,\times \,{\mathrm{CO}}_{2}}^{1/3}=0.17\end{eqnarray} \tag{ 25 }$$

which he based on Equation (12), which in turn can be rewritten as

$$\begin{eqnarray}&&C\,{\Pr }_{\mathrm{air}}^{1/3}=0.068.\end{eqnarray} \tag{ 26 }$$

It is unclear how the factor 0.17 was derived from this, as Pr and Sc of gases are close to one.

Table 2 summarizes the results for the universal prefactor, and in this work we adopt C = 0.14.

Table 2.  Prefactor for Convective Flux from a Horizontal Plate

Source	$C\,{\mathrm{Sc}}^{1/3}$	$C\,{\Pr }^{1/3}$	C
	Mass Transfer	Heat Transfer
Mull & Reiher (1930) a	⋯	(in air)
Jakob (1949)	⋯	0.068 in air
Ingersoll (1970)	0.17 on Mars	⋯
Fujii & Imura (1972) a	⋯	(in water)	0.13
Holman (1990)	Quotes Fujii & Imura (1972)		0.13
Mills (2001)	⋯	⋯	0.14
Hecht (2002)	$0.15\times {0.5}^{1/3}$	⋯	0.15
Incropera et al. (2007)	⋯	x	0.15

Note.

^aOriginal measurements. Ingersoll (1970) is based on Jakob (1949), which fits data from Mull & Reiher (1930).

Download table as:  ASCIITypeset image

Now we consider the values of Sc, D_m, and ν, which are needed for the numerical evaluation of the equations. For an ideal gas $\mathrm{Sc}=\Pr =1$. For a hard-sphere model $\Pr =2/3=0.67$ and $\mathrm{Sc}=5/6=0.83$. Hecht (2002) uses a Schmidt number $\mathrm{Sc}=0.78\times 0.64=0.50$, where 0.78 is the Prandtl number for Mars atmospheric composition and at room temperature. Hence, expression (23) amounts to $\mathrm{Sh}=0.12\,{\mathrm{Gr}}^{1/3}$ for Mars conditions and near melting.

Hudson et al. (2007) includes a comprehensive literature review for values of the vapor diffusivity. Based on this survey, ${D}_{{\rm{H}}2{\rm{O}}.\mathrm{CO}2}=0.15$ cm² s⁻¹ at standard temperature and pressure is a consensus value, which can then be scaled by pressure and temperature. For ideal gases, all three types of diffusivity (κ, D_m = D₁₂, and ν) change with temperature and pressure proportionally to T^3/2/p₀.

The temperature dependence of dynamic viscosity η (Washburn 2003) is

$$\begin{eqnarray}&&\eta ={\eta }_{0}\displaystyle \frac{{T}_{0}+b}{T+b}{\left(\displaystyle \frac{T}{{T}_{0}}\right)}^{3/2}\quad {\rm{(Sutherland}}\,{\rm{formula).}}\end{eqnarray} \tag{ 27 }$$

For CO₂, T₀ = 273 K, η₀ = 13.7 μPa s, and b = 240 K. The dynamic viscosity of a gas is nearly independent of pressure.

With these values for D_m and η, the Schmidt number for CO₂ at 0° C comes out to about 0.47, similar to the value used by Hecht (2002). Both are low compared to theoretical expectation mentioned above. If we adopt a compromise value of $\mathrm{Sc}=0.6$, then

$$\begin{eqnarray}&&C^{\prime} =C\,{\mathrm{Sc}}^{1/3}\approx 0.12\end{eqnarray} \tag{ 28 }$$

for ${{\rm{H}}}_{2}{\rm{O}}$ in CO₂ . This is lower than the 0.17 used by Ingersoll (1970).

Measured sublimation rates in a 7 mbar CO₂ atmosphere (Sears & Moore 2005; Moore & Sears 2006) are consistent with Ingersoll's formula, but those authors did not provide any unitless constants, so a precise comparison of their laboratory measurements with the theoretical prefactor is not immediately available.

Fujii & Imura (1972) measured heat loss from an inclined plane and found that for an inclination angle θ, $\mathrm{Gr}\,\times \,\Pr$ is to be replaced with $\mathrm{Gr}\,\times \,\Pr \,\times \,\cos \theta$. In other words, the convective flux is lower by a factor of ${(\cos \theta )}^{1/3}$ compared to a horizontal plate. For a slope of 45°, this amounts to a factor of ∼0.9, so it is a small effect.

The partial pressure, p_w, can never exceed the total pressure, p₀, so evaporative cooling must compensate any value of heat input. In other words, in the limit ${p}_{w}\to {p}_{0}$, the sublimation rate must diverge, ${E}_{c}\to \infty$, whereas in Equation (18) ${\rm{\Delta }}\rho /\rho \to {M}_{1}/{M}_{w}-1\approx 1.44$ and E_c in (19) remains finite. For p_w ≈ p₀ the analogy between heat and mass transfer breaks down, because it assumes p_w ≪ p₀, so it is justified to change the asymptotic limit of E_c. For p_w ≪ p₀, Equation (18) reduces to $(1-{M}_{w}/{M}_{1}){p}_{w}/{p}_{0}$, so it may as well be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\rho }{\rho }=\displaystyle \frac{{p}_{w}}{{p}_{0}-{p}_{w}}\left(1-\displaystyle \frac{{M}_{w}}{{M}_{1}}\right),\end{eqnarray} \tag{ 29 }$$

which gives almost identical results for small p_w, but diverges at p_w = p₀, as desired.

So far it was assumed that ice sublimates into a dry atmosphere. When sublimation is into an already humid atmosphere, the buoyancy forces are reduced. Near melting this effect is negligible, because, as described above, the humidity of the atmosphere is several orders of magnitude lower than that near the ice. Nevertheless, it is helpful to incorporate this into the equations so the transition to the case of stable ice, when the temperature is below the frost point, becomes continuous. Moreover, no buoyancy-driven convection occurs at low temperature, because wind-caused shear is higher, so the detailed quantitative form does not even matter. Denote the vapor density of the atmosphere as ${\rho }_{\infty }$. When the prefactor ρ_w in Equation (19) is changed to ${\rho }_{w}-{\rho }_{\infty }$, the loss vanishes for ${\rho }_{w}={\rho }_{\infty }$, as desired.

In conclusion, the proposed new parameterization is

$$\begin{eqnarray}\begin{array}{rcl}{E}_{c} & = & 0.14{D}_{m}({\rho }_{w}-{\rho }_{\infty })\\ & & \ \times \ {\left[\displaystyle \frac{g}{{\nu }^{2}}\displaystyle \frac{{p}_{w}}{{p}_{0}-{p}_{w}}\left(1-\displaystyle \frac{{M}_{w}}{{M}_{1}}\right)\cos \theta \right]}^{1/3}{\mathrm{Sc}}^{1/3}.\end{array}\end{eqnarray} \tag{ 30 }$$

It updates the prefactor, diverges at p_w = p₀, and includes the inclination effect.

Figure 2 compares the parameterizations of Ingersoll (1970) and the updated version for a specific atmospheric pressure and a common set of material constants. The offset is due to the adjustment of the prefactor $C^{\prime}$ from 0.17 to 0.12. The most significant difference is the divergence of the updated parameterization when the partial ${{\rm{H}}}_{2}{\rm{O}}$ pressure approaches the total atmospheric pressure.

Zoom In Zoom Out Reset image size

The Grashof scale γ may be calculated as

$$\begin{eqnarray}&&\displaystyle \frac{1}{\gamma }={\left[\displaystyle \frac{g}{{\nu }^{2}}\displaystyle \frac{{p}_{w}}{{p}_{0}-{p}_{w}}\left(1-\displaystyle \frac{{M}_{w}}{{M}_{1}}\right)\right]}^{1/3}.\end{eqnarray} \tag{ 31 }$$

This section provides an overview of the model. Mathematical and implementation details are included in the User Guide that accompanies the online code archive (Schörghofer 2019). (The driving program used in this study is cratersQ_mars_full.f90, release 1.1.6.)

Shadowing by nearby topography (terrain shadowing) defines local horizons. Horizons for each surface element are determined with azimuth rays, with 2° azimuthal resolution, and the highest horizon in each direction is stored. For the purpose of horizon calculations, the topography is represented by triangular surface elements. The horizon elevations and view factors determined from these geometric calculations are stored in a file.

As the Sun moves through the sky, the model then simulates the time evolution of insolation (incoming solar radiation) and surface temperature, using the horizons and view factors as input. The surface energy balance is integrated over time at steps of 1/50 of a solar day (sol) for six Mars years. Surface temperature and insolation are updated at every time step.

The equation governing the energy balance on the surface is

$$\begin{eqnarray}&&{Q}_{\mathrm{total}}+k\displaystyle \frac{\partial T}{\partial z}=\epsilon \sigma {T}^{4},\end{eqnarray} \tag{ 32 }$$

where Q_total is the total absorbed irradiance, k thermal conductivity, T temperature, z depth below surface, emissivity, and σ the Stefan–Boltzmann constant. The equation of heat conduction is solved with a Crank–Nicolson scheme on a grid with spatially varying spacings.

The total radiative energy input is

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{\mathrm{total}} & = & (1-A)({Q}_{\mathrm{direct}}+{Q}_{\mathrm{refl}}+{Q}_{{\rm{a}},\mathrm{sw}}F)\\ & & +\ \epsilon ({Q}_{{\rm{a}},\mathrm{lw}}F+{Q}_{\mathrm{IR}}),\end{array}\end{eqnarray} \tag{ 33 }$$

where A is albedo. The flux Q_direct is determined from the decl. of the Sun, latitude, and hour angle. Due to topography, reflected sunlight (Q_refl) and thermal emission (Q_IR) from other surfaces also need to be added. The terms ${Q}_{{\rm{a}},\mathrm{lw}}$ and ${Q}_{{\rm{a}},\mathrm{sw}}$ are the long-wavelength and short-wavelength contributions from the atmosphere, based on a 0-dimensional model. F is the view factor of the sky.

The surface reemits radiation in all directions (Lambertian), but receives additional heat from surfaces in its field of view. These irradiances depend on view factors V,

$$\begin{eqnarray}&&{Q}_{\mathrm{refl}}=\int \int A^{\prime} \left({Q}_{\mathrm{direct}}^{{\prime} }+{Q}_{\mathrm{refl}}^{{\prime} }\right)V(x,y,x^{\prime} ,y^{\prime} )\,{dx}^{\prime} {dy}^{\prime} ,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{IR}}=\int \int \left[\epsilon \sigma {T}^{{\prime} 4}+(1-\epsilon ){Q}_{\mathrm{IR}}^{{\prime} }\right]V(x,y,x^{\prime} ,y^{\prime} )\,{dx}^{\prime} {dy}^{\prime} ,\end{eqnarray} \tag{ 35 }$$

where primed variables are evaluated at $(x^{\prime} ,\,y^{\prime} )$ and unprimed variables at (x, y). The integrals are over all surface elements with a direct line of view. This is the computationally most expensive component of the model calculations, and limits the investigation to small spatial domains.

The model site is chosen at a latitude of 30° S. At the equator, seasonal effects would be small, and near the poles, peak temperatures would be low, so a mid-latitude location is ideal. A southern latitude is chosen, because the seasons are more extreme in the southern than in the northern hemisphere. Assumed are a surface albedo of 0.12 and an infrared emissivity of 0.98. The chosen albedo value is lower than the global average albedo of Mars, which favors higher peak temperature, but is realistic for many sites. The standard model run assumes a thermal inertia of 400 tiu (the thermal inertia unit, tiu = Jm⁻² K⁻¹ s^−1/2).

The CO₂ frost albedo is 0.65. The CO₂ sublimation temperature is set to 145 K, as appropriate for the southern highlands. The albedo of the water frost is assumed to be the ambient albedo, as the frost is easily darkened by even a thin sublimation lag.

The three-dimensional surface energy calculations are initially carried out without any water frost. At locations of interest, the time series of absorbed energy is stored and the subsurface thermal model is rerun and reequilibrated with this energy as input to the surface energy balance. The presence of water frost is implied by subtracting the latent heat of its sublimation from the energy budget. Far from the melting point, that latent heat (the evaporative cooling) is negligible, but near the melting point it is significant and dependent on the total atmospheric pressure. The geoid-defined zero level of Mars Orbiter Laser Altimeter topography occurs at an average pressure of 520 Pa at L_s = 0° (Smith & Zuber 1998). The 610 Pa isobar lies, on average, 1600 m below that. The lowest-lying basin is the floor of Hellas basin, which lies about 7150 m below the topographic datum. There, the maximum annual pressure reaches 1.2 kPa (Haberle et al. 2001). Below, pressures of 500 Pa and 1 kPa will be considered for the evaluation of the evaporative cooling.

Figure 3 shows model results for a bowl-shaped crater with a depth-to-diameter ratio of 1:5. Annual peak temperatures are above 273 K over the entire domain (on the ice-free surface). A small amount of CO₂ frost forms seasonally on the pole-facing slope. Water frost accumulates continuously for up to hundreds of sols (the temperature continuously remains below the frost point temperature during this period). However, at least 95 sols pass between the end of continuous water frost accumulation and the first time 273 K is reached, so this geometry is not suitable for the melting of frost.

Zoom In Zoom Out Reset image size

For a boulder, here idealized as a half-sphere that sticks out from the surface (Figure 4), the situation is far more favorable. Beyond the southern (poleward) end of the boulder, water frost continuously accumulates for even longer periods, decimeters of CO₂ frost accumulate, and peak temperatures are well above the melting point. Figure 5 shows the time dependence of a location behind the boulder, compared to that of the flat unobstructed surface. The third pixel from the edge was chosen, to avoid any potential corner artifacts. The location is seasonally shadowed around the winter solstice. After the Sun rises again the CO₂ mass soon begins to shrink. After the crocus date, the surface temperature rapidly increases, and within 1 1/4 sols it reaches the melting point—when evaporative cooling is not considered. In other words, the diurnal maximum temperature goes from 145 to 273 K in such a short time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These calculations show that even CO₂ frost can be expected at locations that subsequently reach temperatures above 273 K, consistent with the results by Kossacki & Markiewicz (2004), and this prolongs the duration of ${{\rm{H}}}_{2}{\rm{O}}$ frost. Once the Sun rises, the CO₂ ice begins to sublimate, but the CO₂ –${{\rm{H}}}_{2}{\rm{O}}$ ice composite cannot warm until all of the CO₂ ice has disappeared. By this time, the insolation will be even more intense. After this point, the ${{\rm{H}}}_{2}{\rm{O}}$ frost, none of which has yet been lost to the atmosphere, will experience rapid warming.

The peak surface temperatures behind the boulder are no higher than for an unobstructed surface (Figure 5), so the topography provides no energetic advantage during peak periods. The significant role of the topography is to cast a seasonal shadow. The depression (Figure 3) and the protruding topography (Figure 4) both cause seasonal CO₂ accumulation, but only for the protruding topography (the boulder) is there a rapid transition from cold conditions to low solar incidence angles ("high Sun").

With evaporative cooling, the surface does not reach the melting point (Figure 5). In this model calculation, 1000 Pa and pure ice (no lag deposit) are assumed. The albedo of the ${{\rm{H}}}_{2}{\rm{O}}$ frost is assumed to be only 0.12, as reasonable if darkened by a lag. On Earth, the albedo of old snow is typically 0.4. But a lag deposit only a few particles thick can darken the surface tremendously. On the first full sol after the crocus date, the temperature rises to 256 K and 0.1 kg m⁻² of frost are lost until it first reaches this temperature. This is the equivalent of a 100 μm thick layer of water. The next day, the peak temperature is 260 K and at this point 0.5 kg m⁻² of frost have been cumulatively lost. The evaporative cooling is too strong to allow 273 K to be reached. Melting temperature would be reached, if the loss rate were artificially reduced by a factor of 20.

For an atmospheric pressure of 500 Pa, the values are similar. The peak temperature on the first full sol is 254 K and 0.15 kg m⁻² of frost (0.15 mm thickness) were lost. On the next sol, the peak temperature is 257 K.

The Viking 2 Lander observed almost continuous early frost 10–20 μm thick, and later patchy frost probably 100–200 μm thick (Svitek & Murray 1990). More frost may accumulate in well-shadowed alcoves. Hence, the mass lost within a sol or two after the crocus date is within the amount that can be expected to be present.

Thermal model calculations were carried out not only for a thermal inertia of 400 tiu, but also 200 and 1000 tiu. The different thermal inertias change the crocus date, but the peak temperatures are about the same and melting temperatures were not reached.

A layer of dry material overlying the ice acts as a diffusion barrier and reduces mass loss and evaporative cooling. Such a layer could form as a sublimation lag if dust is included during deposition. The dust content of terrestrial snow is rather low, but deposition rates on Mars are much lower, and therefore the relative dust content could be much higher. If the overlying layer is thin enough, it dampens the diurnal temperature amplitude only slightly so that the underlying ice experiences similar peak temperatures as the surface. Next, we quantitatively evaluate the reduction in sublimation loss due to an overlying dust layer.

When the partial pressure of the water vapor is much lower than the total pressure, the flux through the porous layer is in the form of molecular diffusion,

$$\begin{eqnarray}&&J=-{D}_{{\ell }}\displaystyle \frac{\partial {\rho }_{v}}{\partial z}\approx {D}_{{\ell }}\displaystyle \frac{{\rho }_{w}-{\rho }_{s}}{\zeta },\end{eqnarray} \tag{ 36 }$$

where D_ℓ is the vapor diffusivity of the porous layer. The water vapor density on the surface is ρ_s, while ρ_w is the vapor density at the buried ice surface and equals the saturation vapor density. The approximation, which amounts to a steady-state model, should be valid for a thin layer subjected to nearly periodic temperature cycles.

With advection, there is a pressure difference across the layer, and the flux becomes (Hudson et al. 2007)

$$\begin{eqnarray}&&J=-{D}_{{\ell }}\displaystyle \frac{{\rho }_{0}}{1-c}\displaystyle \frac{\partial c}{\partial z},\end{eqnarray} \tag{ 37 }$$

where c = ρ_v/ρ₀. This relation is no longer useful at c = 1, but it establishes that for high vapor content, advection enhances the flow, and it has the desired property that the flux diverges for ρ_v = ρ₀. The pressure difference can at most be the weight of the overlying layer, which is small, e.g., 1300 kg m⁻³ × 3.7 m s⁻² × 0.01 m = 48 Pa for a 1 cm thick layer. Hence, the partial ${{\rm{H}}}_{2}{\rm{O}}$ pressure still cannot surpass the atmospheric pressure by a significant amount. A practical approximation for the advection correction that retains the qualitative behavior is to replace (36) with

$$\begin{eqnarray}&&J\approx \displaystyle \frac{{D}_{{\ell }}}{1-c}\displaystyle \frac{{\rho }_{w}-{\rho }_{s}}{\zeta }.\end{eqnarray} \tag{ 38 }$$

The flux from the surface to the atmosphere is analogous to (21)

$$\begin{eqnarray}&&J=C^{\prime} {D}_{m}\displaystyle \frac{{\rho }_{s}-{\rho }_{\infty }}{\gamma }.\end{eqnarray} \tag{ 39 }$$

Eliminating ρ_s from (38) to (39) yields

$$\begin{eqnarray}&&J=\displaystyle \frac{{\rho }_{w}-{\rho }_{\infty }}{\tfrac{\zeta }{{D}_{{\ell }}}(1-c)+\tfrac{\gamma }{C^{\prime} {D}_{m}}}=\displaystyle \frac{1}{\tfrac{\zeta }{{D}_{{\ell }}}\tfrac{1-c}{{\rho }_{w}-{\rho }_{\infty }}+\tfrac{1}{{E}_{c}}}.\end{eqnarray} \tag{ 40 }$$

Hence, the diffusion barrier becomes significant when

$$\begin{eqnarray}&&\zeta \gtrsim \displaystyle \frac{{D}_{{\ell }}}{{D}_{m}}\displaystyle \frac{\gamma }{C^{\prime} }\displaystyle \frac{1}{1-c}.\end{eqnarray} \tag{ 41 }$$

When the pore sizes are larger than the mean-free path in the atmosphere, D_ℓ =  D_m, where is the porosity of the dry layer. The mean-free path in the Martian atmosphere is about 10 μm (Hudson et al. 2007). The cross-section for molecular diffusion is restricted by a factor of relative to the entire cross-sectional area available for molecular diffusion in the atmosphere. In this case, the diffusion barrier becomes significant when ζ ≳  γ/C'. The porosity is  ≈ 0.42 for a randomly packed granular medium, and $C^{\prime} \approx 0.12$, so in this case ζ ≳ 3.5γ.

For micron-sized particles (commonly referred to as dust-sized) the pore sizes are smaller than the mean-free path, and D_ℓ/D_m can take on smaller values. For example, for a pore size of 1 μm, ten times less than the mean-free path in the atmosphere, the threshold would be ζ ≳ 0.4γ. Micron-sized particles are known to be suspended in the Martian atmosphere. In conclusion, for the dust layer to be a significant barrier to vapor diffusion, its thickness needs to be on the order of the Grashof length scale γ.

Figure 6 shows values for γ according to Equation (31) with a linear vertical axis. At 260 K, γ is about one centimeter. For an overlying layer to significantly reduce the mass flux, and therefore, evaporative cooling, it would have to be nearly that thick. Dust is incorporated in frost and snow, but to form a sublimation lag, a far thicker layer of ice must first be lost. If the dust content is 1%, then about 10 cm of ice need to be lost to leave behind a 1 mm thick layer of dust. A protective dust layer could also originate from a dust storm.

Zoom In Zoom Out Reset image size

Seasonal CO₂ frost is easily decimeters thick. If enough dust is incorporated into the CO₂ ice during accumulation, a dust layer will build up as the CO₂ sublimates back into the atmosphere, leaving behind water frost and the dust.

A typical diurnal thermal skin depth for dust on Mars is 4 cm, so a layer a few mm thick will not significantly dampen the temperature amplitude.

The thermal model calculations in Section 3 demonstrated that sudden transitions from frost-accumulating conditions to near-melting conditions occur, but ultimately there is not enough energy available to compensate for the evaporative cooling. Here, we focus on the energetics, with and without a protective layer, using an idealized energy input.

An idealized form for the absorbed solar flux is

$$\begin{eqnarray}&&{Q}_{\mathrm{total}}=0.96\times (1-A)\displaystyle \frac{{S}_{0}}{{R}^{2}}\cos i,\end{eqnarray} \tag{ 42 }$$

where A is albedo, S₀ = 1360 W m⁻² the solar constant, R the distance from the Sun in au, and i is the incidence angle, assumed to change linearly with time. The factor of 0.96 represents net absorption in the Martian atmosphere. This approximation represents the Sun rising and setting at an equatorial location at equinox. The equilibrium surface temperature T is calculated by numerically solving the nonlinear equation ${Q}_{\mathrm{total}}=\epsilon \sigma {T}^{4}+{L}_{{\rm{H}}2{\rm{O}}}J(T)$, where L_H2O is the latent heat of sublimation and J the sublimation rate (40).

Figure 7 shows the surface temperature and the time-integrated loss over a quarter of a sol (from morning until noon). At an atmospheric pressure of 500 Pa, below the triple point, 273 K cannot be reached. For exposed ice, the surface temperature reaches 264 K and about 1.2 kg m⁻² ice, if available, is lost during this quarter sol. With 2 mm of very fine sublimation lag, these values remain almost the same (267 K and 1.1 kg m⁻²). For an atmospheric pressure of 1 kPa and without overlying layer, the peak temperature is 269 K and the loss 1.1 kg m⁻² by noon. For an atmospheric pressure of 1 kPa and a 3 mm thick layer of micron-sized particles, 273 K is reached and the ice loss until this point of time is about 0.5 kg m⁻² (about 500 μm). This last result is for rather favorable conditions: solar energy input at the equator, at perihelion, equilibrium temperature (zero thermal inertia), a high pressure of 1000 Pa, and a thick lag consisting of micron-sized particles. Realistically, this ideal combination cannot be expected for any place on Mars, so the conclusion is that crocus melting of pure ice does not occur on present-day Mars.

Zoom In Zoom Out Reset image size

For a frost albedo of 0.4 instead of 0.15, the noontime temperatures are lower by about 5 K in all four cases. Having enough refractory material to darken the surface contributes notably to the peak temperature.

These results show: (1) the consequences of evaporative cooling are severe. For all but the most favorable conditions, the melting point is not reached. (2) A sublimation lag of a few mm reduces the sublimation loss only slightly. At 500 Pa, a 2 mm sublimation lag of micron-size particles increased the peak temperature only by 3 K. (3) Peak temperatures within about 10 K of the melting point within two sols of the crocus date are realistic. The full thermal model calculations with exposed (but darkened) ice resulted in peak temperatures of 254–260 K within the first two sols after the crocus date. For the energetically more favorable situation considered in this section, exposed ice reaches 264–269 K at noon. For realistic energetics and a thin sublimation lag, 263 K is a realistic peak temperature. (4) Since the energetics considered in this section is close to optimal for the present orbit (at perihelion, Sun rising to zenith, dark frost), changes in latitude and axis tilt will not suffice to cause crocus melting. The three-dimensional thermal model calculations demonstrated that a sudden transition from frost-accumulating conditions to hot conditions can occur. They were carried out for only one site, but other site locations or site parameters would still not have a higher maximum energy input than was assumed for the idealized model.