Stability of mid-latitude snowpacks on Mars

Icarus 196 (2008) 565–577

Contents lists available at ScienceDirect

Icarus www.elsevier.com/locate/icarus

Stability of mid-latitude snowpacks on Mars K.E. Williams a,∗ , O.B. Toon a , J.L. Heldmann b , C. McKay b , M.T. Mellon a a b

Laboratory for Atmospheric and Space Physics (LASP UCB 392), Department of Atmospheric and Oceanic Sciences, University of Colorado, Boulder, CO 80309-0392, USA NASA Ames Research Center, Division of Space Sciences and Astrobiology, Moffett Field, CA 94035, USA

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 30 November 2006 Revised 25 January 2008 Available online 8 April 2008

Christensen [2003. Nature 422, 45–48] suggested that runoff from melting snowpacks on martian slopes might be responsible for carving gullies. He also suggested that snowpacks currently exist on Mars, for example on the walls of Dao Valles (approximately 33◦ S). Such snowpacks were presumably formed during the last obliquity cycle, which occurred about 70,000 years ago. In this paper we investigate a specific scenario under conditions we believe are favorable for snowpack melting. We model the rate at which a snowpack located at 33◦ S on a poleward-facing slope sublimates and melts on Mars, as well as the temperature profile within the snowpack. Our model includes the energy and mass balance of a snowpack experiencing diurnal variations in insolation. Our results indicate that a dirty snowpack would quickly sublimate and melt under current martian climate conditions. For example a 1 m thick dusty snowpack of moderate density (550 kg/m3 ) and albedo (0.39) would sublimate in less than two seasons, producing a small amount of meltwater runoff. Similarly, a cleaner snowpack (albedo 0.53) would disappear in less than 9 seasons. These results suggest that the putative snowpack almost certainly could not have survived for 70,000 years. For most of the parameter settings snowpack interior temperatures at this latitude and slope do reach the melting point. Under most conditions melting occurs when the snowpack is less than 10 cm thick. The modeled snowpack will not melt if it is covered by a 1 cm dust lag. In general, these findings raise interesting possibilities regarding gully formation, but perhaps mostly during a past climate regime when snowfall was expected to have occurred. If there currently are exposed snowpacks on martian mid-latitude slopes, then these ice sheets cannot last long. Hence they might be time variable features on Mars and should be searched for. © 2008 Elsevier Inc. All rights reserved.

Keywords: Ices Mars, surface

1. Introduction Gully-like features have been noted on Mars since the Mars Orbiting Camera (MOC) onboard the Mars Global Surveyor (MGS) arrived in 1997 (Malin and Edgett, 2000). Gullies typically are about 1 km long (Heldmann and Mellon, 2004). Many of the gully-like features appear to be geologically young (within the last few million years) since most are not cratered and many are superimposed on surfaces which themselves are geologically young. Several models have been proposed for the formation of gullylike features (hereafter just referred to as gullies). Some researchers (Musselwhite et al., 2001) propose that a liquid subsurface CO2 reservoir may be responsible for carving the gullies. There are numerous problems with the CO2 flow hypothesis, however, including the fact that the flow-channel morphologies are inconsistent with the physics of CO2 flow (Heldmann and Mellon, 2004). Other researchers (Mellon and Phillips, 2001) have argued that subsurface aquifers with impermeable boundaries may be re-

*

Corresponding author. E-mail address: [email protected] (K.E. Williams).

0019-1035/$ – see front matter doi:10.1016/j.icarus.2008.03.017

©

2008 Elsevier Inc. All rights reserved.

leased by fractures in the ice plug, causing a cascade of liquid that could have incised some gullies. Heldmann (Heldmann et al., 2007; Heldmann and Mellon, 2004) studied hundreds of gullies and determined that most of the gullies were best explained by the characteristics of subsurface fluid sources rather than surface snowmelt, ground ice melt or wind. Following the work of previous researchers (Clow, 1987; Farmer, 1976) Christensen (2003) suggested that gullies formed when exposed dusty snowpacks slowly melt. The meltwater becomes runoff and incises gullies underneath the melting snowpack. Specifically, he hypothesized the following: 1. Water is transported from the martian poles to mid-latitudes in the form of snow during periods of high obliquity. 2. Melting occurs at mid-latitude during low obliquity producing liquid water that is stable beneath snow. 3. Meltwater causes surface erosion and the formation of small gullies. 4. Gullies incised into the substrate are observed where snow has been removed. 5. Patches of snow remain today where they are protected from sublimation by a layer of desiccated dust/sediment.

566

K.E. Williams et al. / Icarus 196 (2008) 565–577

Fig. 1. Proposed snowpack in Dao Valles. This photo is a portion of MOC image M09-02875 and covers an area approximately 1.5 × 3.0 km located at approximately 33.3◦ S, 92.9◦ E (Christensen, 2003). The arrow points to a proposed snowpack, along with a possible exposed gully. North is towards the top.

Fig. 2. Eccentricity and obliquity for the last 90,000 years for Mars. These data were supplied by Laskar (private communication). Note that the current obliquity (25.19◦ ) was last achieved approximately 70,000 years ago when the eccentricity was at a temporary maximum.

6. Melting could be occurring at present in favorable locations in these snowpacks. Christensen noted images of a curious mantling in the Dao Valles region (∼33◦ S by 93◦ E) as evidence (Christensen, 2003). He suggested that this mantling (Fig. 1) is a remnant snowpack/icecemented soil. He pointed out that several gullies can be seen protruding from this mantling, indicating that they therefore could be the result of the snowpack melting. There are numerous unresolved issues with the above hypotheses. First, are the mantled deposits in Fig. 1 actually dusty snowpacks? If so then they must have survived for a very long time, since according to Christensen (2003) they were presumably laid down when the obliquity was substantially higher than at present, probably more than 70,000 years BP (Fig. 2). Second, can melting occur below the surface of mid-latitude snowpacks, and if so can it do so under current conditions and in enough volume to create runoff and subsequent erosion? Hypothesis #2 (above) as-

serts that melting occurs at mid-latitude and that the liquid water would be stable beneath the snow. There is an important difference, however, between liquid water being produced and liquid water being produced in sufficient quantities to achieve runoff and erosion (Clow, 1987; Paterson, 1994). For water to achieve runoff capability there needs to be a refrozen (impermeable) snow boundary along the base (which may even inhibit erosion). The presence of such a layer has been modeled by Clow but his results suggest that runoff was achieved only under a very narrow set of circumstances: runoff was achieved only if the snowpack was thin (∼23 cm) and very dusty (1000 ppmw) at 7 mbar ambient pressure (Clow, 1987). In other situations the liquid water may be briefly stable but never achieve runoff capability before it is refrozen in the basal layer of the snowpack. Clow (1987) found that if the snowpack was too thin (less than ∼23 cm) or too thick (greater than ∼35 cm) then runoff was not produced. If the snowpack was too thin the meltwater was used entirely in creating the refrozen basal snow layer and no runoff was produced. If the snowpack was too thick there was not enough internal heating to produce meltwater, since the heat is more efficiently wicked away from warm spots in thicker snowpacks. In addition, according to Clow (1987), reliable and widespread runoff and melting would need substantially higher atmospheric pressure than currently exist on Mars. It is now doubtful whether these higher atmospheric pressures ever existed, since recent research suggests that the permanent polar caps (a CO2 reservoir) would themselves probably only contribute 0.36 mb of additional CO2 if they completely sublimated (Byrne and Ingersoll, 2003). More generally, however, Clow’s investigation made several simplifications that need to be addressed. Clow (1987) did not explicitly include mass loss in his model (so his snowpack never disappeared). He also did not explicitly treat the sensible heat flux, latent heat flux and atmospheric heat flux terms in the surface energy balance but instead parameterized the sum of these terms to be sinusoidally varying around a mean value (derived from midlatitude values). He did not apply his model to different latitudes and slopes/aspects, but rather to mid-latitudes on level ground. Lastly, his radiative transfer model did not allow for inhomogeneous dust mixtures within the snowpack. Christensen (2003) relies heavily on Clow’s calculations and it is necessary to examine the applicability of Clow’s work to slopes. Note that in this study we are specifically modeling snow, and not bubble-free ice (so-called “blue” or “black” ice). The subsurface melting possibilities of bubble-free ice are intriguing, but the likelihood of bubble-free ice on martian slopes is unknown, and it seems more likely that ice or snow on martian would contain both bubbles/pore spaces as well as dust. In this paper we examine the question of whether a 1 m thick snowpack located in the southern mid-latitudes and on a slope such as in Fig. 1 could melt under current Mars climate conditions. We also investigate whether a 1 m thick snowpack formed at high obliquity on the same slope could survive to the present day and whether such a snowpack could melt and/or provide runoff (which could form a gully). We focus on conditions exemplified by Dao Valles gullies where apparently 1–10 m of mantle deposit exists (Christensen, 2003). We use a model similar to the model of Clow (1987) but with applicability to slopes and with atmospheric terms developed in much greater detail than in Clow’s model, as well as a more sophisticated radiative transfer component. 2. Description of our model The full model solves heat, radiation and mass transfer equations in the ice. Below we first discuss the heat transfer equation within the snow. We then describe the boundary conditions on the heat balance equation. Following that we describe the mass bal-

Snowmelt on Mars

Table 1 Parameter

Low value

Base case

High value

Snow density Atmospheric pressure Surface atmos. precipitable water pressure Wind speed Snow dust content Latitude Slope inclination Slope aspect Orbital eccentricity Obliquity Ground thermal conductivity (W/(m K)) Ground volumetric C p (J/(m3 K)) Emissivity of snow Areothermal heat flux (mW/m2 )

400 kg/m3 7 mbar 0.037 Pa

550 kg/m3 8.75 mbar 0.056 Pa

700 kg/m3 11 mbar 0.074 Pa

0 m/s 10 ppmw 30.0◦ S 0◦ 0◦ 0.05 15◦ 0.045

2.5 m/s 100 ppmw 33.0◦ S 21◦ 135◦ 0.09 25.19◦ 0.172

5 m/s 1000 ppmw 42.0◦ 40◦ 270◦ 0.13 35◦ 0.60

1.297 0.95 30

1.381 0.97 30

1.507 0.99 30

567

integrate Eq. (1) and obtain the following for a layer k of thickness  z:

 ∂ Φ =− ( J k+1/2 − J k−1/2 ) + Ψ  z, ∂t

where Φ has been normalized per unit area and Ψ is a source term for the volume element. Since Eq. (2) is for a control volume, the node k corresponds to the center of the kth volume element, and k − 1/2, k + 1/2 correspond to the upper and lower faces of the volume element. In the case of our 1-D snow model, the conserved quantity is the layer energy U normalized per unit area where U = ρ C p T z. The energy U of the layer is then evolved in time by solving:

     j +1/2 ∂T  ∂T  ∂U j −1/2

+ F net − F net = − k( T )  − k( T )  ∂t ∂ z j +1/2 ∂ z j −1/2 + ζ  z.

Fig. 3. The energy terms used in our model.

ance equation. Finally we discuss the radiative transfer equations. Table 1 summarizes the parameter values we have chosen, and the ranges of variation we considered in sensitivity tests. The model is iterated over time until the remaining snowpack is approximately 2 cm thick. When only ∼2 cm of snow remains we stop the model and declare the snowpack to be obliterated. In order to initialize the snowpack temperature profile, the annual average snowpack surface temperature was computed for the slope scenario. The snowpack and soil substrate was then initialized with the resulting temperature. Various other initialization scenarios were tried, and this scenario was by far the most computationally efficient. 2.1. Energy transfer within the snow layer Fig. 3 shows the energy terms of the model. Our model uses a control volume approach, where for a given volume V we can equate the time rate of change of the product of density ρ (in kg/m3 ) and a conserved quantity φi with a surface integral of a flux averaged over a specified time interval. In general:

∂ ∂t



  φi ρ dV = − J i · dS + ψ dV , i



S

(2)

(1)

V

where V ψ dV is a generic source term. Here i is a general index, S denotes the surface of the volume element V , and J a flux density. Given the usual assumption of homogeneity of the conserved quantity over the designated control volume, we may

(3)

Here the overbar indicates that the fluxes are averaged over an interval greater than or equal to the designated integration timestep interval, and since in our case U is normalized per unit area, it has units of J/m2 , giving the left-hand side of Eq. (3) units of W/m2 . The last term on the right-hand side of Eq. (3) (ζ  z) is the average energy source/sink rate contributed by the arrival/departure of percolated meltwater in the layer for the current timestep, and has units of W/m2 . F net is the net solar flux at the control volume boundary due to solar radiation penetrating the snow as described in Section 2.4. The temperature-dependent coefficient of thermal conductivity for snow, k( T ), is computed using the Brailsford and Major model (Clow, 1987). The Brailsford and Major model combines thermal conductivity across snow pore spaces with the thermal conductivity along ice grains. It formulates an effective thermal conductivity which is comprised of the thermal conductivity of CO2 gas across the pores, the thermal conductivity associated with water-vapor diffusion across the pores, the radiative conductivity across pore cavities and the thermal conductivity of solid ice along the ice grains. The heat capacity of the snow (C p ) was computed using a temperature-dependent form for ice (Paterson, 1994), C p ( T ) = 152.5 + 7.122T





in J/(kg K) .

(4)

The snow consists of coarse granular snow, often referred to as firn. The snow density, ρ , for the base case (Table 1) was chosen to be 550 kg/m3 , since this density was surmised by Clow to be the result of several seasons worth of densification (metamorphosis) (Clow, 1987). Fluffy and high density snows were considered in sensitivity tests. A density of 400 kg/m3 corresponds to freshlyfallen fluffy snow (Clow, 1987). A high snow density of 700 kg/m3 corresponds to very wet terrestrial snow (Paterson, 1994). The vertical computational grid used in our model is a variable layer mass grid where layer thickness is diagnosed based on densities of the layer constituents. Most layers are 1 cm thick but thicknesses are allowed to vary throughout the snowpack as mass is gained or lost in various layers in the snow column. Clow used a logarithmic grid to allow for very thin layers at the surface (∼1 mm) and thick layers in the interior, but we elected to keep ∼1 cm layers since the penetration of the sunlight is generally much larger than 1 cm (Clow, 1987). We also chose 1 cm layers because we wanted to resolve <10 cm snowpacks adequately, while still being computationally efficient. The snowpack rests on a soil substrate (i.e. the ground). The ground is also subdivided into layers whose thickness is 1 cm. Note that not only do we consider the ground beneath the ice as just discussed, but we also solve an energy balance equation for the ice free ground near the ice sheet, since that ground is likely to determine the local air temperature.

568

K.E. Williams et al. / Icarus 196 (2008) 565–577

Ground thermal conductivity was varied between 0.045 and 0.60 W/(m K) with a base case value of 0.172 W/(m K), which is in the range of values for dry mineral soils (Usowicz et al., 2006). Clow uses a value of 0.30 W/(m K), where he found that most of his melting occurred. He considers this to be a low value for conductivity and found that low ground conductivity favors melting. In order to use a very conservative value (favorable for melting) we therefore chose the very low ground thermal conductivity of 0.172 W/(m K). As explained later, however, we did not find much model sensitivity to the ground conductivity. Ground volumetric C p varied between 1.297 to 1.507 MJ/(m3 K) with a base case value of 1.381 MJ/(m3 K). Clow used 1.28 MJ/(m3 K), which is somewhat lower than our value. Note, however, that our results are largely insensitive to variations in ground C p . 2.2. Energy boundary conditions The substrate boundary condition for the snow layer energy equation (3) is a constant upward areothermal heat flux of 30 mW/m2 , though estimates vary by a few mW/m2 (Schubert et al., 1992). The sublimation of the snowpack is not sensitive to this tiny heat flux, so we did not consider a range of areothermal heat fluxes in our sensitivity tests. The surface energy balance in our model consists of the following components (Fig. 3): infrared emission from the snow, solar absorption by the snow, atmospheric infrared emission to the snow, infrared emission from an adjacent surface to the snow, sensible heat flux and latent heat flux due to mass loss. The energy for the surface layer is evolved by the following expression:

  ∂U ∂T  = −k( T )  − εσ T 4 + ( F −1/2 − F +1/2 ) ∂t ∂ z bottom  ∂M , + A H + S H + F IR − L ∂t

(5)

where it is understood that the right-hand side terms are each averaged over a time interval at least as large as the integration timestep. The choices of various parameters in this equation are summarized in Table 1, and discussed below. The quantity U = mC p T is the area-normalized energy of the surface layer of ice/dust/liquid mixture of mass M and has units of J/m2 , hence the left-hand side of Eq. (5) has units of W/m2 . Again the equations are solved on a mass grid; when required, thickness  z is simply diagnosed from the layer mass since the relative mass contributions of liquid water, dust and ice are tracked for each layer. The F IR term is the infrared heat flux density contribution of a planar surface at the foot of the slope, and is described later in this article. The infrared emission term, εσ T 4 includes the emissivity ε of the snow surface. Some terrestrial models set the emissivity of all permanent snow features to be 0.95 (Box et al., 2004). The NCAR Land Surface Model uses a value of 0.97 for snow surface emissivity (Bonan, 1996). Others used a value of 0.99 for Mars (Heldmann and Mellon, 2004). The values for snow IR emissivity in sensitivity tests vary between 0.95 and 0.99, and we chose the base case to be the average of these two values (0.97). Since the mass and energy in a layer is evolved in time, the temperature which is diagnosed for a snow layer corresponds to the center of a layer. Hence the “surface” temperature in our model is actually about 5 mm below the snow “surface,” depending on the current thickness of the surface layer. The concept of the snow “surface” is problematic, however, given that the radii of the modeled snow particles themselves are 4 mm, and their packing arrangement at the surface is unknown. Indeed the surface may have additional roughness features other than the snow particles themselves, which when combined with the particles yields an

even larger overall roughness scale of the order of centimeters (or more). We believe, therefore, the depth of the surface temperature measurement to be a negligible issue since our modeled snow particles and surface roughness uncertainty are comparable to the depth variation of the temperature measurement. For the surface layer (Eq. (5)), the flux term k( T ) ∂∂Tz is evaluated at the bottom boundary of the surface layer. The F −1/2 − F +1/2 term is the average net solar flux absorbed by the surface layer (direct and scattered radiation) for the timestep. The flux is computed for a transparent atmosphere overlying the snow. The solar flux is a function of time, latitude and slope inclination/aspect. It has been computed here using fixed astronomical variables (fixing obliquity and eccentricity, and using present-day longitude of perihelion), since the timescales for the sublimation turn out to be relatively short. Note that using present-day longitude of perihelion is favorable for southern-hemisphere melting. A latitude of 33◦ S was assumed corresponding to the approximate location of Dao Valles (Fig. 1). Recent findings have suggested that gullies are commonly found between 30◦ S and 42◦ S, hence the sensitivity tests varied latitude between these latitudes (Heldmann et al., 2007). We have estimated a slope aspect of 135◦ , and slope inclination angle of 21◦ ± 4◦ . Slope aspect is measured clockwise from north, and slope inclination angle is measured up from the horizontal plane. The slope inclination angle was varied from 0◦ to 40◦ in order to capture the range between “level” ground and the angle of repose. This slope inclination range basically covers the range of alcove slopes reported by Heldmann and Mellon (2004). In order to assess model sensitivity to aspect, we varied the aspect between 0◦ and 360◦ . Note that 0◦ aspect means a north-facing slope, a 90◦ aspect means eastward-facing slope and 180◦ aspect means a due-south facing slope. We have used the equivalent slope concept (Dingman, 2002) to account for slope inclination and aspect. This concept accurately represents the solar flux on a plane parallel semi-infinite slope. It makes no approximations other than ignoring shadowing and other three-dimensional issues. The technique is essentially a latitude and longitude shift; a slope of 30◦ poleward is equivalent to a 30◦ poleward shift in latitude. East/west slopes become a longitude adjustment of their sunrise/sunset times, as well as their solar noon times. For example, a due-east slope would have the same sunrise time as a flat surface, but its sunset is earlier than on a flat surface. The new sunset time is the same as that of a flat surface located further east at a different longitude (hence the concept of an equivalent longitude). The equivalent latitude and the adjusted longitude are:





ϕ = sin−1 cos( I )∗ sin( L 0 ) + sin( I )∗ cos( L 0 )∗ cos( A ) , Ωadj =

tan−1 (sin( I )∗ sin( A )) cos( I )∗ cos( L 0 ) − sin( I )∗ sin( L 0 ) ∗ cos( A )

,

(6) (7)

where ϕ is the equivalent latitude, I is the slope inclination, L 0 is the current latitude and A is the slope aspect. Omega is the equivalent longitude. The new sunrise and sunset times for the slope are then computed by: sunrise_t = sunset_t =

− cos−1 (− tan(ϕ )∗ tan(δ)) − Ωadj

ω − cos−1 (− tan(δ)∗ tan( L 0 ))

ω

,

,

(8) (9)

where δ is the declination and ω is the angular velocity of Mars. The solar noon and current solar time is computed by: solar_noon =

−Ωadj

ω

,

solar_t = local_time − (12.33 + solar_noon),

(10) (11)

Snowmelt on Mars

Fig. 4. Average daily insolation over a flat surface for several latitudes in the northern hemisphere of Mars. These curves are for current obliquity and eccentricity values.

569

Fig. 5. Average daily insolation over a flat surface for several latitudes in the southern hemisphere of Mars for current obliquity and eccentricity values.

where 12.33 corresponds to the time in (terrestrial) hours of local noon for Mars. The hour angle is then: h = ω∗ solar_t .

(12)

The cosine of the solar incidence angle, i, on the slope is then computed by: cos _i = cos(ϕ )∗ cos(δ)∗ cos(h + Ωadj ) + sin(δ)∗ sin(ϕ ).

(13)

F d is finally computed as: F d (t ) = S 0 (r0 /r )2 cos _i ,

(14)

where r is the current Sun–Mars distance, r0 is the mean Sun– Mars distance (over a year). An annual average value of 588.98 W/ m2 (Kieffer et al., 1992) was used for S 0 . Orbital eccentricity for the base case is set at the current value of 0.0934, and varied between 0.05 and 0.13 since these values were achieved within the past 500 my (Laskar, J., 2004. Personal communication: Martian obliquity and eccentricity values for 500,000 years BP). Obliquity is varied between 15◦ and 35◦ for similar reasons, with the base case using the current value of 25.19◦ . In our model we have ignored the visible wavelength optical depth of the atmosphere. In essence, ignoring the dust opacity in our model will maximize the possibility of melting by maximizing insolation. Figs. 4 and 5 show average daily insolation over a level ground at several latitudes, using current eccentricity and obliquity values. As can be seen on the figures, high latitudes receive no insolation for a considerable part of a martian year, but they receive a great deal of insolation (more than 200 W/m2 ) during the summer. Mid-latitude insolation over flat topography during the year varies between 50 and 150 W/m2 in the winter to close to 200 W/m2 in the summer. Fig. 6 shows the daily average insolation for various polewardfacing slope angles at a latitude of 33◦ S. As can be seen in the graph, slope angle has a noticeable impact on the amount of insolation reaching the surface. Varying the slope angle from 0◦ to 30◦ can reduce the wintertime insolation at 33◦ S latitude by almost 50 W/m2 per day. Fig. 7 shows the daily average insolation for various slope aspects. Note that an aspect of 180◦ corresponds to a poleward-facing slope. During the wintertime, it can be seen in the graph that aspect can also drastically affect the daily-averaged amount of insolation incident on our area of interest. East-facing slopes will experience the local noon (halfway between local sunrise and local sunset) earlier than a flat slope would experience

Fig. 6. Average daily insolation over a variety of slope inclination angles (slopes are facing 135◦ , which is east of poleward) at a latitude of 33◦ S for Mars. This latitude and aspect corresponds approximately to the northern wall of Dao Valles, the location of the hypothesized snowpack deposits.

Fig. 7. Average daily insolation over several slope aspects (slopes all have inclination 21.0◦ ) at a latitude of 33◦ S for Mars. This latitude and slope inclination corresponds approximately to the northern wall of Dao Valles, the location of the hypothesized snowpack deposits. Note that 90◦ aspect corresponds to due east and 180◦ aspect is poleward facing.

570

K.E. Williams et al. / Icarus 196 (2008) 565–577

local noon because local sunset occurs earlier. West-facing slopes experience local sunrise shifted later in the day. We should reiterate that we are considering a flat slope, with a flat snowpack on the surface of the slope. We envision that actual extensive snow fields, wider than the gully at top, are needed to supply enough liquid to be of interest in carving the gully. We are not considering the melting possibilities of snowpacks inside of gullies. The 3-D geometry of gullies are perhaps worthy of study, but there are an infinite number of gully geometries. In addition, a flat slope with a snowpack on the surface will provide more meltwater (to incise a gully) than a small snowpack in the bottom of a gully. Also, we do not consider shadowing effects, given the infinite variety of scenarios, which make it extremely difficult to treat adequately. Also we believe many gullies are on slopes with essentially unobstructed views of the sky, where shadowing effects are minimal. If a crater has gullies on a rim, then the height of a neighboring rim subtends only a tiny fraction of the visible sky from the perspective of the gully location: if a 1 km deep crater is 10 km wide, and the gully is located halfway down the crater rim, then the opposite crater wall occupies <7% of the sky (∼(arcsin(1/10))/90 = 5.7/90 < 7%). Nevertheless, we have included a potential source of radiation which our slope may “see,” namely surface re-radiation and reflection, albeit in a limited form. Consider an infinite plane at the foot of the model slope, stretching out to a flat horizon. Our modeled slope will intersect a small amount of reflected shortwave solar radiation from the plane, as well as a small amount of radiated longwave IR. Accordingly, we have used the appropriate equations from Mellon and Phillips (2001). They include the absorbed solar shortwave radiation reflected from the plane at the foot of the slope:





F vis = (1 − A α ) A 0 F s (0.5) 1 − cos(α )

(15)

and the absorbed infrared radiation which is emitted by the plane at the foot of the slope:





F IR = εα ε0 σ T 04 (0.5) 1 − cos(α ) ,

(16)

where (0.5)(1 − cos(α )) accounts for the mean incidence angle of the plane onto the slope. Here F s is the incident solar flux on a flat surface, ε is the infrared emissivity, T is the surface temperature and A the albedo. The slope surface and plane surface are indicated by α and 0 subscripts, respectively. The F terms (in Eqs. (3) and (5)) are the net solar fluxes. We compute these explicitly using a radiative transfer numerical model, rather than specifying an ice albedo as discussed below. There are two reasons for this choice. First we wished to make the albedo a function of the dust content of the ice. Secondly we wished to include solar heating that occurs with depth in the ice. The downwelling atmospheric infrared flux term A H is parameterized in our model to be: A H CO2 = 3.0 × 105 ln(1 + 1.97 × 10−3 P )

× exp(−960.0/ T s ) ergs/cm2 /s,

(17)

where P is the ambient atmospheric pressure, T s is the surface temperature, and σ is the Stefan–Boltzmann constant. The parameterization is from Toon et al. (1980) and Gierasch and Goody (1968) but scaled by 0.5. We scaled by 0.5 in order to fit our atmospheric heating to the amount implied by the atmospheric radiative transfer model and temperature profile of McKay and Davis (1991), which was derived from the work of Pollack et al. (1987). In Fig. 8 we show our fit to the McKay and Davis data. We obtain a value for the downward flux of approximately 7.2 W/m2 for a CO2 atmosphere, a surface temperature of 200 K and an ambient pressure of 875 Pa.

Fig. 8. Here we show the calculated average martian surface temperature using our atmospheric heating parameterization (dashed line) and a fit by McKay and Davis (1991) to a 1-d radiative transfer model result (solid line). We differ by ∼1 K over the pressure range of our model runs, which is acceptable for our purposes.

Atmospheric pressure P was set to 875 Pa since Viking Lander 2 at 48◦ N latitude found a yearly average of 875 Pa with variation between 700 and 1000 Pa (Zurek, 1992). The sensitivity tests considered a pressure range between 700 and 1100 Pa. The sensible heat flux term SH is the sum of a free-convection (buoyancy) term and a forced (wind-driven) term. The forced term is: SH1 = ρ C p A d u ( T a − T s ),

(18)

where A d is a drag coefficient, u is the mean wind speed, T a is the temperature of the atmosphere just above the surface and T s is the surface (ice) temperature. The drag coefficient is poorly constrained, and we have chosen the value of 0.002 as recommended in Paterson (1994) as the value for the drag coefficient over an icy surface. We have elected not to vary this parameter in a sensitivity test because it is multiplied by wind speed, and wind speed is included in our sensitivity tests where it is varied between 0 and 5 m/s. Viking lander surface wind speeds at 25◦ N latitude suggest the wind speeds at that location are often between 0 and 5 m/s (Murphy et al., 1990). The free-convection sensible heat flux term (Holman, 2002) is: SH2 = 0.15( T a − T s )k(Gr Pr)1/3 ,

(19)

where Gr and Pr are the Grashof and Prandtl numbers, respectively (defined in Appendix A), and k is the thermal conductivity for CO2 gas (Holman, 2002). The sensible heat flux term SH is then SH = SH1 + SH2 . The surface air temperature T a (a necessary parameter in many of the above equations) is computed by solving the energy balance model over dry sandy soil and taking the air temperature to be the ground surface temperature. This seems reasonable given that the snowpacks are not extensive, and the air temperature over the snow is most likely controlled by air advected there from the surrounding soil-covered areas. The final term in the surface energy equation (Eq. (5)) is the latent heat times the mass loss. L is the latent heat of sublimation of ice, found by adding the latent heats of fusion and vaporization. The mass loss term is discussed in the next section. 2.3. Mass loss equations We considered snow-column mass loss to occur both at the surface of the ice (as sublimation) and at the base of the snow column (as runoff). The model surface layer is unique compared to

Snowmelt on Mars

other layers in that it is permitted to lose mass both by sublimation and by percolation. We modeled percolation of water through the snowpack by a simple threshold rule, where the meltwater in a layer was moved down to the next lower level (in the next timestep) if the liquid mass fraction of meltwater was greater than or equal to 5% for the layer. This threshold was also used by Clow (1987). The percolated meltwater is allowed to refreeze if/when it encounters a sufficiently cold snow layer. If the meltwater reaches the base (substrate), it is allowed to leave the snow column and is considered to be runoff. No soil substrate infiltration is modeled. Also, because of the 1-D nature of the model, no lateral movement of meltwater is permitted. The sublimation mass loss process (the mass loss mechanism from the snow surface) was modeled as the sum of two terms; following the work of Hecht (2002), the sublimation mass balance term ∂∂Mt is comprised of two terms:

∂M = m1 + m2 . ∂t Following the work of others (Ivanov and Muhleman, 2000) we use the following equation for forced convection mass loss as: m1 = − M w

1 kT

Au (e − e sat ),

(20)

where M w is the molecular weight of water, k is the Boltzmann constant, T is the snow surface temperature, A is a drag coefficient for icy surfaces which is taken to be 0.002 in our model as recommended in Paterson (1994), e is the partial pressure of water vapor, and e sat is the saturation vapor pressure of water (Ivanov and Muhleman, 2000; Paterson, 1994). We compute the free convection mass loss term using: m2 = 0.17(e − e sat ) M w

1 kT

 D

ρ

ρ



g

ν2

1/3 ,

(21)

where

ρ

ρ

= 26e sat /(44P − 26e sat ),

where ν is the kinematic viscosity of CO2 gas, e sat is the water saturation vapor pressure, R is the gas constant and D = 32T 3/2 / P is the diffusion coefficient of H2 O in CO2 (Ingersoll, 1970). The derivation of Eq. (21) can be found in Appendix A. The parameterization for e sat depends on the temperature. For temperatures above 273 K we used the Lowe and Ficke parameterization (Seinfeld and Pandis, 1998). For temperatures between 170 and 273 K we used the parameterization of Marti and Mauersberger (1993). For temperature less than 170 K we used the parameterization of Mauersberger and Krankowsky (2003). A reasonable parameterization for e (in mbar) for water is difficult to formulate given the interannual and annual variability in the column water vapor abundance (Haberle and Jakosky, 1992). However, we have taken the surface water abundance shown in General Circulation Model (GCM) simulations for 33◦ S and 25◦ obliquity (Mischna et al., 2003). According to Mischna et al. (2003), Mars GCM simulations under current obliquity suggest that at southern mid-latitudes, zonally averaged water vapor varies throughout the year between 10 and 20 pr μm (precipitable microns). We have therefore parameterized water vapor pressure as the pressure corresponding to a constant mixing ratio, giving a column 15 pr μm of water (annual average). We did not include seasonal changes in the water abundance, even though Mischna et al. (2003) found some variation. There is probably considerable interannual and annual variability in the column water vapor abundance (Haberle and Jakosky, 1992). Therefore for sensitivity testing we varied this value between 10 and 20 pr μm corresponding to the range that Mischna et al. (2003) found in their GCM simulations for current Mars orbital conditions.

571

Laboratory data (Hecht, 2002) have shown the Ingersoll (1970) equation (Eq. (21)) to be correct in the no-wind situation, however it is not clear from the literature precisely how to combine the free and forced convective mass loss equations in a low wind situation (in our case, a 2.5 m/s mean wind speed). However, we believe his equation is appropriate in our wind regime, given the following general argument: Assume a 100 m wide snow patch, and a uniform friction velocity u ∗ of 0.5 m/s. Then the wind at the surface is able to cross the patch in 200 s. Assume a mass loss rate similar to that of Fig. 9, or ∼1.0 kg/day. Assume an average temperature of the snow surface to be 260 K. Using the diffusion coefficient of water in CO2 of 153 cm2 /s, and a moisture roughness length for the snowy surface of 10 cm we find a 10 cm thick layer of air over the snow would become saturated in less than 1 s. Since the layer of the surface becomes saturated so quickly relative to the time required to blow an air parcel over the snow patch, we believe the buoyancy effects to still be in effect even though the wind is blowing, and therefore Ingersoll’s equation is relevant to the modeled snow patch. Fig. 9 shows sublimation rates and surface temperatures for a very dusty snowpack during a summer day on a 21◦ slope (and 135◦ aspect) at 33◦ S lat. Note that the temperature and sublimation maxima occur at approximately 12:00 local time even though the slope is facing toward the morning Sun. As shown in Fig. 10, insolation on a southeast facing slope peaks before noon, however the temperature response of the snow is slightly delayed. Free convective mass loss is the dominant mass loss process. 2.4. Radiative transfer equations Our model uses the shortwave radiative transfer code of McKay et al. (1994) which is a two-stream code that is based on the algorithms in Toon et al. (1989). Our radiative transfer code fully models a non-homogeneous dust distribution within the snow. Most current models of martian ice sheets specify the ice albedo and assume that the light does not penetrate beyond their first model level in the ice. Our goal is to make the snow albedo consistent with the dust content, and also to determine if the light penetration matters to the snow sublimation rate. Clow (1987) modeled the snow albedo to be consistent with the dust content of the snow, and also suggested that light penetration might be important in distributing the energy deep within the snow. Clow (1987) developed a radiation code based on analytic solutions to the two stream equations for this purpose. Since it is an analytic solution the dust is assumed to be homogeneously distributed throughout the snow layers. We have used the code of McKay et al. (1994) because homogeneity may not be a good assumption, since the top of the snow will become dust-enriched as the snow sublimates. However, as dust becomes concentrated in the top layer it is expected that solar radiation will eventually be completely absorbed in the top model layer. It is possible that radiative transfer in layers below the top layer is no longer important when dust concentrations become high. Nevertheless we chose to model the full non-homogeneous snowpack for completeness sake. As in Clow (1987), the dust content of the snow was varied between 10 and 1000 parts per million by weight (ppmw). The base case value we used was 100 ppmw, which corresponds to slightly dusty snow with an albedo of 0.39. This value of snow albedo was very close to that used by Mischna et al. (2003) as an annual average value for the North Polar ice cap albedo for their simulations (0.37). An upper value of 1000 ppmw was specifically chosen because Clow (1987) found that maximal melting (and, more importantly, runoff ) occurred at this dust content. Completely clean snow is not expected to occur on Mars. A layer of pure dust was computed to have a visible wavelength albedo of 0.13. The average value of the surface albedo in

572

K.E. Williams et al. / Icarus 196 (2008) 565–577

Fig. 9. Sublimation rates and snow surface temperature during a typical Mars spring day ( L s = 180◦ ) on a 21◦ slope at 33◦ S latitude. Slope aspect is 135◦ . This figure is for a dusty snowpack with 100 ppmw of dust.

Fig. 10. Insolation for a typical spring day ( L s = 180◦ ) for our slope, located at 33◦ S, with slope angle 21◦ and slope aspect 135◦ . The insolation peaks on our slope around 10 a.m., but as can be seen in Fig. 9 there is a delayed temperature response.

the region considered in this paper is ∼0.15 (Putzig, N.E., 2006. Personal communication: Thermal inertia and albedo estimate), so our value of 0.13 is conservative since it maximizes the chance of melting dusty snow by maximizing the air temperature. The high dust value of 1000 ppmw produces a computed snow albedo of 0.25, which corresponds to very dusty snow. The dust properties used in our numerical model follows Clow (1987). The dust is red with a fixed 2 μm grain size. The dust grains were assumed to be spherical. Both the single scattering albedo parameterization and asymmetry factor parameterization for the dust are discussed in Clow (1987). In addition, the radiative transfer numerical model includes reflection from the soil substrate at the base of the snow column.

A sublimation lag on the surface would block a significant amount of sunlight as the surface layer becomes dust enriched, effectively shutting off sunlight penetration at depth. Our results suggest that ∼95% of the downward solar flux is attenuated by the top 1 cm of dusty snow if the dust content at the surface layer is 10,000 ppmw. If a 1 m snowpack starts with 1000 ppmw dust, then by the time 4 cm of snowpack are left the surface layer has enriched itself to approximately 95,000 ppmw. However, for snowpacks that are not very dusty the light can penetrate deeply. For example with the base case dust, 16% of the downward solar flux reaches a depth of 10 cm. We have included a lag deposit in our numerical model by monitoring the dust in the topmost layer of snow that is freed by ice sublimation.

Snowmelt on Mars

Table 2 Sensitivity test of our numerical model for a 1 m thick snowpack Parameter that was varied

Resulting snowpack lifetime (years, days)

Meltwater runoff (l/m2 )

Highest T reached (K)

Snow density 400–700 (kg/m3 ) Atmos. pressure 7–11 (mbar) Atmos. precip. water 0.037–0.074 (Pa) Wind speed 0–5 (m/s) Snow dust content 10–1000 (ppmw) Latitude 30 to 42 (◦ S) Slope inclination 0, 10, 30, 40 (◦ )

1 y, 377 d–2 y, 314 d 1 y, 356 d–1 y, 380 d 1 y, 370 d–1 y, 370 d

0.73–1.0 0.8–1.3 1.4–1.4

273.15–273.15 273.15–273.15 273.15–273.15

2 y, 267 d–1 y, 325 d 8 y, 425 d–1 y, 283 d

0.2–1.3 0.0–0.9

273.15–273.15 269.7–273.15

1.6–0.0 0.4, 0.1, 0.0, 0.0 0.0, 1.9, 1.4, 4.5 0.0–0.0 0.0–9.1 2.2–0.3

273.15–273.15 273.15, 273.15, 273.15, 267.8 273.15, 273.15, 273.15, 273.15 273.15–272.17 272.7–273.15 273.15–273.15

y, 368 d–1 y, 371 d

0.9–1.0

273.15–273.15

y, 339 d–1 y, 423 d

2.5–0.0

273.15–270.8

1 0 4 Slope aspect 0–90–180–270 (◦ ) 0 1 Eccentricity 0.05–0.13 3 Obliquity 15–35 (deg) 5 Ground thermal conductivity 1 0.045–0.60 (W/(m K)) Ground volumetric C p 1 1.297–1.507 (J/m3 K) Snow emissivity 0.95–0.99 1

y, y, y, y, y, y, y, y,

341 d–2 y, 358 d 421 d; 1 y, 219 d; 290 d; 16 y, 350 d 485 d; 1 y, 328 d; 288 d; 1 y, 320 d 288 d–1 y, 210 d 292 d–0 y, 361 d 371 d–1 y, 365 d

3. Numerical model results In this section we present the results of our numerical model. We first examine the base case, where we have run the model using parameters initialized with values that are the most reasonable (justified in Section 2). We then consider sensitivity to the various parameters by varying the uncertain parameters between reasonable limits. 3.1. Results for base case A base case was run under current climate conditions using fixed parameter values which were justified in the model description section above (Table 1). We also provide results in Table 2 for a range of eccentricity and obliquity conditions, though a complete and thorough examination of paleo-conditions remains an open research question. Fig. 11 shows the time to lose 1 m of snow is less than two seasons (starting at L s = 90◦ ). The maximum (including the surface layer) temperature of the snowpack reaches the melting point (273.15 K) for the base case, but only when the snowpack is 5.9 cm thick. The base case does produce a small amount of meltwater runoff (∼1.4 l/m2 ). The model was run at 10-min timesteps in order to capture diurnally-varying quantities without resorting to diurnal averages. Unless stated otherwise, seasonally and diurnally varying quantities either were computed explicitly or as 10-min (model timestep) averages, rather than using diurnal or seasonal averages. As outlined previously, the model was initialized by setting the temperature to be isothermal with the annual average snowpack surface temperature. For all of the cases reported herein, the model simulation was started at L s = 90◦ , which corresponds to southern winter solstice (presumably the time at which the southern hemisphere snow deposits occurred). Fig. 12 shows the temperature profiles for the hottest time of day for several days during the first martian year. As can be seen in Fig. 11, the snowpack considered here lasts less than two seasons before completely sublimating and melting. This snowpack begins with a relatively high albedo of 0.39. The snow does reach melting temperatures high enough for melting to occur. The melting for the base case (1 m) occurs during summertime and when the snowpack is relatively thin (5.9 cm). At that time, the surface layer dust content has been considerably enriched to approximately 5800 ppmw, which causes a high heating rate to

573

occur in the surface layer. At the same time, the ice is sufficiently warm due to both the favorable season (summertime) and the fact that the substrate has been allowed to warm (due to the favorable season and sunlight absorption). The melting occurs at the surface layer, and spreads downward as warm water percolates to the underlying layers. A thicker snowpack (10 m) melts in a similar manner when it gets very thin. However, the thicker snowpack also melts and produces small amounts of runoff in summertime when it is approximately 3.7 m thick (or less), and its surface layer has enriched to approximately 54,000 ppmw. At that time the underlying ice is also warm enough for melting to occur at the surface layer. From this we may conclude that for melting to occur there needs to be both a high heating rate at the surface as well as relatively warm ice underneath. The interior ice is warm when the season is favorable and/or when the substrate is warm. We have also determined the minimum constant amount of water that must be present in the atmospheric column in order to keep an exposed 1 m snowpack, of the type specified in the base case, stable. For current Mars orbital conditions we find that the atmospheric water pressure must be ∼7.5 Pa (implying a column of ∼200 pr μm) for an exposed clean snowpack (snow density 550 kg/m3 , 100 ppmw dust, albedo 0.39) to remain stable under current climate conditions. 4. Results for sensitivity tests Table 2 lists the parameters that were varied in order to test the sensitivity of the model. The sensitivity tests were done for a 1 m thick snowpack. The model snow lifetime was most sensitive to slope inclination, followed by dust content of snowpack, obliquity and eccentricity of the orbit and density of the firn. From Table 2 it appears that the snowpack melting and runoff amounts are also sensitive to various model parameters. However, upon further analysis, our model findings suggest that the amount of meltwater runoff is not an important indicator of snowpack lifetime. Instead, the production of meltwater is closely tied to whether the snowpack was sufficiently thin during the late spring and early summer to produce melting. If the snowpack happens to be too thick during a particular summer season, no melting or runoff will be produced. Such timings are more or less fortuitous, and have little to do with the actual lifetime of the snowpack. As can be seen in Table 2, the snowpack lifetime appears to be sensitive to snow density. For a thin (1 m) snowpack of low density (400 kg/m3 ) the lifetime was less than two seasons. For higher density snow (700 kg/m3 ) the lifetime is less than three seasons. While the lifetime difference is largely attributed to the greater mass in the snow column for the high density case, the difference is also due to the variation in thermal conductivity of snows of various densities. Snow density plays a role in the thermal conductivity of the snow, and the thermal conductivity of the snow affects how quickly heat is carried away from hot spots within the snowpack. A lower bulk density of the snow results in a lower thermal conductivity. A lower bulk density for snow of 400 kg/m3 may be possible from snow fall. However, the rates of snow metamorphosis are almost certain to be large for martian atmospheric conditions. Clow (1987) estimated that within one or two seasons of the snowfall, metamorphosis would convert the (initially) fluffy snow to very granular snow of ice grain radii of 4 mm and snow bulk density of ∼550 kg/m3 . If, as Christensen (2003) hypothesized, the snow fell during a previous high obliquity cycle, the snowfall presumably fell at least 70,000 years BP. It is likely that our base case estimate of snow density is conservative (from a melting standpoint), and that in fact a snowpack that was wellpreserved would be even more dense than we estimate for our base case.

574

K.E. Williams et al. / Icarus 196 (2008) 565–577

Fig. 11. The base case: Daily snow thickness as a function of time (Mars years) for latitude of 33◦ S, slope of 21◦ and 135◦ aspect under current Mars orbital conditions. The simulation was started at L s = 90◦ . The snow has a density of 550 kg/m3 and 100 ppmw dust content. The entire 1 m snowpack disappears in less than 2 seasons. The max daily surface temperature data have been smoothed.

Fig. 12. Snowpack temperature profile for three different days during the first year of model run. The location is the base case location of 33◦ S, 21◦ slope and 135◦ aspect. All three profiles are from noon local time.

The lifetime of the snowpack depends on the amount of dust in the snow. A greater amount of dust within the snowpack decreases the bolometric albedo. Also, greater dust loading provides more opportunity for solar absorption by the dust grains and thereby increases the temperature (at depth) of the snowpack. The effect of depositing energy at depth in the snowpack is to increase the temperature, especially if the thermal conductivity of the snowpack is low. The mass loss sensitivity to slope angle, latitude, obliquity and eccentricity follow from the effects of these parameters on the insolation. As can be seen in Fig. 13, lower obliquities result in lower daily average summertime insolation. It can be seen in Figs. 6 and 7 that insolation can vary considerably depending on aspect

and slope angle. It is interesting that a very steep slope (40◦ ) can extend the snowpack lifetime to approximately 17 seasons, given the severely reduced insolation on such a slope. Eccentricity affects insolation as well. As can be seen in Fig. 14, higher eccentricity results in higher daily average insolation during the summer season. 5. Discussion and conclusion Given the snow model we have used, it is not possible that a 1 m thick, dirty snowpack persisted at 33◦ S and 21◦ slope with 135◦ aspect for longer than two seasons. We conclude that the mantled deposit in Fig. 1 is not likely to be an exposed snowfield since it would have a very short lifetime. Melting temper-

Snowmelt on Mars

Fig. 13. Daily averaged insolation for various obliquities over flat ground, at 33◦ S latitude.

Fig. 14. Daily average insolation for three different eccentricities over flat ground. Latitude is 33◦ S. Longitude of perihelion was 250.87◦ and longitude of vernal equinox was 109.29◦ .

atures are reached when the snowpack was very thin (5.9 cm), however, and a small amount of meltwater runoff was produced (∼1.4 l/m2 ). Note that this slope geometry is specific to the northern wall of Dao Valles. There may have been other times in the past which favored snow buildup simply because there may have been a great deal of water in the atmosphere. For example, the obliquity was close to 33◦ approximately 500,000 years BP, and there must have been much more atmospheric water at that time. Recent gcm simulations suggests that 500 pr μm of atmospheric water were present when the obliquity was 45◦ (Mischna and Richardson, 2005). As mentioned previously, such atmospheric water amounts would probably keep the snowpack stable. Nevertheless, as our model indicates, it is very doubtful that an exposed snowpack could have survived for more than a few seasons at midlatitudes unless water vapor concentrations were extremely large. There is a scenario where the snowpack could survive: immediate burial by a thick dust lag. We were able to completely arrest the melting of a 1 m snowpack by placing a 1 cm dust lag on the snow surface. While melting was not occurring in that scenario, the snowpack still sublimated in less than two seasons. Other modeling has shown that on a similar slope to the one we studied the snowpack could be stable if buried by as little as 10 cm of dust, assuming an annual average of 10 pr μm of water in the atmosphere (Farmer and Doms, 1979). Growing a lag thickness of

575

10 cm is quite attainable by very thick snowpacks (10 s of meters) containing a large amount of dust (1000 ppmw). The thermal skin depth for the soil in our model is ∼5 cm, which suggests that even a 5 cm lag would be enough to protect the underlying snow from the Sun. Such a burial would have to occur very quickly (within a few years of snow deposition) in order for the snowpack to be preserved. However, if lag deposition were to occur quickly (for example, via a slope failure) this would be expected to be a highly localized phenomenon, and probably not extensive enough to blanket large tracts of snowdrifts. In addition, for melting to occur there must also be a rapid exposure of the preserved snowpack. If the snowpack is exposed too gradually, for example by the lag being blown off the snowpack surface, then the snowpack would rapidly sublimate from beneath the lag. Of course under conditions in which snowfalls are stable against sublimation (due to high atmospheric water amounts for instance) a thick snowfield could grow, and then be capped by a dust layer as it sublimated under later conditions. Such a layer could not melt however. In addition to the rapid sublimation of the snowpack we also find that melting and runoff occurs for many of the cases, but mainly when the snowpack is <10 cm thick. While an interesting finding in itself, we believe this raises serious doubts about Christensen’s hypotheses concerning present-day snowpack-fed gully erosion given the fact that the putative snowpacks would have then disappeared long ago. Our findings concerning melting agree fairly well with those of Clow (1987). Clow (1987) found melting and runoff occurred at 7 mbar atmospheric pressure, when the snow was 23 cm thick and contained 1000 ppmw of dust for level surfaces at mid-latitude insolation values. We were able to replicate this result for a level surface, though our melting and runoff occurred when the snowpack was somewhat thinner (∼5.8 cm). Unfortunately, Clow (1987) does not provide enough information for us to precisely replicate his simulations. Consequently we have been unable to precisely compare his findings with our own model results, although they are clearly similar. Modeling by Mellon and Jakosky (1993) has shown that the ice table depth on a slope like the one we have modeled could be as shallow as 10 cm, suggesting that a 10 cm lag could theoretically completely arrest the sublimation of the snowpack. Recent laboratory experiments suggest that a 50 mm lag can slow the sublimation of water ice on Mars by a factor of 10 (Chevrier et al., 2007), but this diffusive loss is still rapid. The depth of the ice table is maintained mostly by the dust layer’s thermal insulation and not by the dust layer’s diffusive insulation (Mellon and Phillips, 2001). Burying ice preserves the ice mostly because the dust keeps the ice colder. It is possible that the addition of salts could result in even greater melting. However, on closer examination the possibility of brines or sulfate salts being mixed into the dust seems unlikely to significantly depress the freezing point, given that we do not have much dust mixed into the snow. For example, the amount of magnesium sulfate required to depress the freezing point of water by 3.6◦ C is 16% (anhydrous solute weight percent) or 32.76% (hydrated solute weight percent) (Weast, 1985). Our average dust concentrations are between 100 ppmw and 1000 ppmw. Even if all of the dust were comprised of salt, the salt concentration would only be 0.01% (100 ppmw salt) or 0.1% (1000 ppmw salt), which is not even enough salt to depress the freezing point by 0.1 K (Weast, 1985). A similar argument can be made regarding sodium chloride, where a 4 K depression would require 6% (anhydrous solute weight percent). Nevertheless it is possible that the lag and/or substrate beneath the snowpack could have lots of salts, thereby depressing the freezing point locally within the snowpack. In conclusion, we find that it is unlikely that mid-latitude gullies are forming today via the mechanism that Christensen posits. Primarily the issue is the lack of a current source of snow. It would

576

K.E. Williams et al. / Icarus 196 (2008) 565–577

require the source of the snow fields be able to deposit about 10 cm of snow each year so that the snow could be melted in the spring or summer. It is not plausible that ancient snow fields (created in a time of greater snow deposition rates, or more prolonged snow deposition) have persisted into modern times to melt now. However, it may be possible that the gullies formed in the past via snowmelt. The loss of snow under current climate conditions occurs much more rapidly than climate change. If any ancient mid-latitude snowpacks exist today they must be stabilized by a dust/soil lag, from which melting is not plausible. Clow’s model was cited by Christensen as providing theoretical evidence for the plausibility of snowmelt forming under current climate conditions. Our model suggests that mid-latitude snowpacks both sublimate and melt very quickly under current Mars conditions. Note, however, that if we assume the premise that partially exposed midlatitude snowpacks currently exist on Mars, then it is likely that melting and runoff is occurring as suggested by Christensen (2003) based on Clow (1987). There are alternative explanations for ice sheets appearing on Mars slopes, such as ice sheets on slopes being formed by springs [as Heldmann and Mellon (2004) and Heldmann et al. (2005) suggested]. Regardless of the specific origin of the ice (should it exist), our model suggests such ice sheets cannot last long. Hence they would be time variable features on Mars and should be searched for.

The authors wish to thank two anonymous reviewers for their helpful comments. K.E.W. also wishes to acknowledge financial support from a NASA GSRP fellowship, and to thank Gary Clow, Mark Williams and W. Tad Pfeffer for helpful discussions. Appendix A. Derivation of equation for free convection mass loss In this appendix we derive the equation for free-convection mass loss. Following the work of Hecht (2002) and Holman (2002) we begin with the basic heat transfer equation between a warm wall and a cold atmosphere: A

= h( T w − T ∞ ),

(A.1)

where q is the heat flux, A is the area (hereafter taken to be 1 m2 ), h is the heat transfer coefficient, T w is the temperature of the wall and T ∞ is the temperature of the atmosphere. The value of h is not easy to determine, and instead we introduce the ratio of convection heat transfer to conduction heat transfer: Nu =

qconvection qconduction

=

h T k xT

=

hx k

,

q = h T = (Nu)k T .

k ρC p

=

μC p k

,

ν2

(A.5)

.

In this case x3 is a unit volume. The parameter β is a volume expansion coefficient, and is defined as:

β=

1 ∂V V ∂T

= =

1 ρ

for mass convection,

ρ T 1 T

(A.6)

for thermal convection,

where in this case the derivative is evaluated at a constant pressure. V is a volume element and T is the gas temperature. The Grashof number is often written as: Gr = x3

 g ( ρ )  ρ . 2

(A.7)

ν

The product of the Prandtl and Grashof numbers is called the Rayleigh number: Ra = Gr Pr.

(A.8)

Nu = 0.15(Gr Pr)1/3 ,

(A.9)

which is appropriate for 8 × 10 < Ra < 1 × 10 . Analogous to Eq. (A.3), mass transfer can be written as: 6

11

∂m C =D Nu, ∂t x

(A.10)

where again x is a length scale taken to be 1 m. For evaporation (mass loss) the Nusselt number is defined as: Nu = 0.15(Gr Le Pr)1/3 ,

(A.11)

where Le is the Lewis number (the ratio of thermal diffusivity to mechanical diffusivity). Using the evaporation form of the Nusselt number in Eq. (A.11), we then have:



E = 0.15(Gr Pr Le)1/3

D C



(A.12)

.

x

Ingersoll (1970) simplifies Eq. (A.12) by approximating: 0.15(Pr Le)1/3 ≈ 0.17, which then gives the final form of the evaporation formula:

 E = m2 = 0.17D C (Gr )1/3 = 0.17D C

(A.3)

Next we introduce three dimensionless numbers that frequently arise when studying fluids. The Prandtl number is defined as the ratio of momentum diffusivity (kinematic viscosity) to thermal diffusivity:

ν

g β T x3

(A.2)

where x is a length scale and k is conductive heat transfer coefficient (which is easier to measure). Take x to be 1 m. We can then write:

Pr =

Gr =

According to Hecht (2002), the Nusselt number for free convection can be written as:

Acknowledgments

q

relative differences between the velocity boundary layer and the thermal boundary layer. Hence when Pr is tiny, heat diffuses relatively quickly. The Grashof number is the ratio of buoyancy forces to viscous forces. It is defined as:

(A.4)

where μ is the dynamic viscosity and k is the thermal conductivity of the material. The Prandtl number is a measure of the

= 0.17(e − e sat ) M w

1 kT

 D

ρ

ρ

g ρ

ν2ρ



g

ν2

1/3 1/3 .

(A.13)

References Bonan, G., 1996. A Land Surface Model (LSM Version 1.0) for Ecological, Hydrological, and Atmospheric Studies: Technical Description and User’S Guide. NCAR, p. 150. Box, J.E., Bromwich, D.H., Bai, L.S., 2004. Greenland ice sheet surface mass balance 1991–2000: Application of Polar MM5 mesoscale model and in situ data. J. Geophys. Res. Atmospheres 109, doi:10.1029/2003JD004451. D16105. Byrne, S., Ingersoll, A.P., 2003. A sublimation model for martian south polar ice features. Science 299, 1051–1053. Chevrier, V., Sears, D.W.G., Chittenden, J.D., Roe, L.A., Ulrich, R., Bryson, K., Billingsley, L., Hanley, J., 2007. Sublimation rate of ice under simulated Mars conditions and the effect of layers of mock regolith JSC Mars-1. Geophys. Res. Lett. 34, doi:10.1029/2006GL028401. L02203.

Snowmelt on Mars

Christensen, P.R., 2003. Formation of recent martian gullies through melting of extensive water-rich snow deposits. Nature 422, 45–48. Clow, G.D., 1987. Generation of liquid water on Mars through the melting of a dusty snowpack. Icarus 72, 95–127. Dingman, S.L., 2002. Physical Hydrology. Prentice Hall, Upper Saddle River, NJ. Farmer, C.B., 1976. Liquid water on Mars. Icarus 28, 279–289. Farmer, C.B., Doms, P.E., 1979. Global seasonal-variation of water-vapor on Mars and the implications for permafrost. J. Geophys. Res. 84, 2881–2888. Gierasch, P., Goody, R., 1968. A study of thermal and dynamical structure of martian lower atmosphere. Planet. Space Sci. 16, 615–646. Haberle, R.M., Jakosky, B.M., 1992. The seasonal behavior of water on Mars. In: Kieffer, H., Jakosky, B.M., Snyder, C.W., Matthews, M.S. (Eds.), Mars. Univ. of Arizona Press, Tucson, pp. 969–1016. Hecht, M.H., 2002. Metastability of liquid water on Mars. Icarus 156, 373–386. Heldmann, J.L., Mellon, M.T., 2004. Observations of martian gullies and constraints on potential formation mechanisms. Icarus 168, 285–304. Heldmann, J.L., Toon, O.B., Pollard, W.H., Mellon, M.T., Pitlick, J., McKay, C.P., Andersen, D.T., 2005. Formation of martian gullies by the action of liquid water flowing under current martian environmental conditions. J. Geophys. Res. Planets 110, doi:10.1029/2004JE002261. E05004. Heldmann, J.L., Carlsson, E., Johansson, H., Mellon, M.T., Toon, O.B., 2007. Observations of martian gullies and constraints on potential formation mechanisms. II. The northern hemisphere. Icarus 188, 324–344. Holman, J.P., 2002. Heat Transfer. McGraw–Hill, New York. Ingersoll, A.P., 1970. Mars—Occurrence of liquid water. Science 168, 972–973. Ivanov, A.B., Muhleman, D.O., 2000. The role of sublimation for the formation of the northern ice cap: Results from the Mars Orbiter Laser Altimeter. Icarus 144, 436–448. Kieffer, H.H., Jakosky, B.M., Snyder, C.W., Matthews, M.S., 1992. The planet Mars: From antiquity to present. In: Kieffer, H., Jakosky, B.M., Snyder, C.W., Matthews, M.S. (Eds.), Mars. Univ. of Arizona Press, Tucson, pp. 1–33. Malin, M.C., Edgett, K.S., 2000. Evidence for recent groundwater seepage and surface runoff on Mars. Science 288, 2330–2335. Marti, J., Mauersberger, K., 1993. A survey and new measurements of ice vaporpressure at temperatures between 170 and 250 K. Geophys. Res. Lett. 20, 363– 366. Mauersberger, K., Krankowsky, D., 2003. Vapor pressure above ice at temperatures below 170 K. Geophys. Res. Lett. 30 (3), doi:10.1029/2002GL016183. 1121. McKay, C.P., Davis, W.L., 1991. Duration of liquid water habitats on early Mars. Icarus 90, 214–221.

577

McKay, C.P., Clow, G.D., Andersen, D.T., Wharton Jr., R.A., 1994. Light transmission and reflection in perennially ice-covered Lake Hoare, Antarctica. J. Geophys. Res. Oceans 99, 20427–20444. Mellon, M.T., Jakosky, B.M., 1993. Geographic variations in the thermal and diffusive stability of ground ice on Mars. J. Geophys. Res. Planets 98, 3345–3364. Mellon, M.T., Phillips, R.J., 2001. Recent gullies on Mars and the source of liquid water. J. Geophys. Res. Planets 106, 23165–23179. Mischna, M.A., Richardson, M.I., 2005. A reanalysis of water abundances in the martian atmosphere at high obliquity. Geophys. Res. Lett. 32, doi:10.1029/ 2004GL021865. L03201. Mischna, M.A., Richardson, M.I., Wilson, R.J., McCleese, D.J., 2003. On the orbital forcing of martian water and CO2 cycles: A general circulation model study with simplified volatile schemes. J. Geophys. Res. Planets 108 (E6), doi:10.1029/2003JE002051. 5062. Murphy, J.R., Leovy, C.B., Tillman, J.E., 1990. Observations of martian surface winds at the Viking Lander-1 site. J. Geophys. Res. Solid Earth Planets 95, 14555–14576. Musselwhite, D.S., Swindle, T.D., Lunine, J.I., 2001. Liquid CO2 breakout and the formation of recent small gullies on Mars. Geophys. Res. Lett. 28, 1283–1285. Paterson, W., 1994. The Physics of Glaciers. Butterworth Heinemann, Oxford. Pollack, J.B., Kasting, J.F., Richardson, S.M., Poliakoff, K., 1987. The case for a wet, warm climate on early Mars. Icarus 71, 203–224. Schubert, G., Solomon, S.C., Turcotte, D.L., Drake, M.J., Sleep, N.H., 1992. Origin and thermal evolution of Mars. In: Kieffer, H., Jakosky, B.M., Snyder, C.W., Matthews, M.S. (Eds.), Mars. Univ. of Arizona Press, Tucson, pp. 147–183. Seinfeld, J.H., Pandis, S.N., 1998. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. Wiley, New York. Toon, O.B., McKay, C.P., Ackerman, T.P., 1989. Rapid calculation of radiative heating rates and dissociation rates in inhomogeneous multiple-scattering atmospheres. J. Geophys. Res. Atmospheres 94, 16287–16301. Toon, O.B., Pollack, J.B., Ward, W., Burns, J.A., Bilski, K., 1980. The astronomical theory of climatic-change on Mars. Icarus 44, 552–607. Usowicz, B., Lipiec, J., Marczewski, W., Ferrero, A., 2006. Thermal conductivity modeling of terrestrial soil media—A comparative study. Planet. Space Sci. 54, 1086– 1095. Weast, R.C., 1985. CRC Handbook of Chemistry and Physics. CRC Press, Boca Raton, FL. Zurek, R., 1992. Comparative aspects of the climate of Mars: An introduction to the current atmosphere. In: Kieffer, H., Jakosky, B.M., Snyder, C.W., Matthews, M.S. (Eds.), Mars. Univ. of Arizona Press, Tucson, pp. 799–817.