Icarus 205 (2010) 269-282 Contents lists available at ScienceDirect ELSEVIER Icarus journal homepage: www.elsevier.com/locate/icarus Evaluating the meaning of "layer" in the martian north polar layered deposits and the impact on the climate connection Kathryn E. Fishbaugh *, Shane Byrne, Kenneth E. Herkenhoff, Randolph L Kirk, Corey Fortezzo, Patrick S. Russell, Alfred McEwen Smithsonian National Air and Space Museum, Center for Earth and Planetary Studies, P.O. Box 37012, MRC 315, Washington, DC 20013, USA ARTICLE INFO Article history: Received 18 August 2008 Revised 14 February 2009 Accepted 4 April 2009 Available online 10 May 2009 Keywords: Mars, Polar Geology Mars, Climate Mars, Polar Caps ABSTRACT Using data from the High Resolution Imaging Science Experiment (HiRISE) aboard the Mars Reconnais- sance Orbiter, we reassess the methods by which layers within the north polar layered deposits (NPLD) can be delineated and their thicknesses measured. Apparent brightness and morphology alone are insuf- ficient for this task; high resolution topographic data are necessary. From these analyses, we find that the visible appearance of layers depends to a large degree on the distribution of younger, mantling deposits (which in turn is partially influenced by inherent layer properties) and on the shape and location of the particular outcrop. This younger mantle partially obscures layer morphology and brightness and is likely a cause of the gradational contacts between individual layers at this scale. High resolution images reveal that there are several layers similar in appearance to the well-known marker bed discovered by Malin, M., Edgett, [(., 2001. J. Geophys. Res. 106, 23429-23570. The morphology, thicknesses (4-8 ? \/2m), and separation distances (5-32 ? V2 m) of these marker beds, as gleaned from a high resolution stereo dig- ital elevation model, lend insight into the connection between stratigraphy and climate. Published by Elsevier Inc. 1. Introduction and background Since their discovery in Mariner 9 (Murray et al., 1973) and Vik- ing images (Cutts et al., 1976), layers within both the north (NPLD) and south polar layered deposits (SPLD) have been characterized broadly as alternating bright and dark stripes (Fig. 1). Many have linked this "striping" to changing relative amounts of dust and ice in the atmosphere as a result of changing orbital parameters and, hence, changing climate (e.g., Murray et al., 1973; Cutts and Lewis, 1982; Laskar et al., 2002; Levrard et al., 2007). But analysis of Mariner 9 images of the SPLD show that layer topography affects the retention of frost and therefore the brightness of individual lay- ers (Herkenhoff and Murray, 1990). Here, we use high resolution image and topographic data from the High Resolution Imaging Science Experiment (HiRISE; (McEwen et al., 2007)) on board the Mars Reconnaissance Orbiter to supplement the preliminary findings of Herkenhoff et al. (2007) by taking a new look at how a "layer" can be delineated, by examining small-scale layer characteristics, and by discussing some of the implications of these analyses for global climate change. Since the PLD likely contain the most complete record of relatively recent climate change on Mars, it is crucially important to obtain a realistic understanding of what a "layer" is and what * Corresponding author. Fax: +1 202 786 2566. E-mail address: fishbaughke@si.edu (K.E. Fishbaugh). 0019-1035/$ - see front matter Published by Elsevier Inc. doi:l 0.1016/j.icarus.2009.04.011 its characteristics and stratigraphic position tell us about the con- temporary martian climate. Delineating layers is a necessary first step toward this understanding. 2. Datasets For all analyses in this study, we use data from HiRISE, including all summertime images of the NPLD acquired during 2006-2007. The HiRISE camera is capable of taking three-color (red: 550-850 nm; blue-green: 400-600 nm; infrared: 800-1000 nm) images; the pixel scale of red images is typically 25-33 cm while the BG and IR images are typically binned 2x2 or 4 x 4. The ~255 x 320 km orbit results in south polar images that are typi- cally 25-27 cm/pixel, while those over the north pole are usually 31-33 cm/pixel. We have also constructed a Digital Elevation Model (DEM) of a trough in the north polar cap, centered near 87.1 ?N, 93.0?E, by ste- reo analysis of HiRISE images PSP_001739_2670 and PSP_001871_2670 (Fig. 2a). Because these images were taken on separate orbits whose ground tracks are not parallel at this high latitude, the overlap area is approximately hexagonal with dimen- sions of 6 x 16 km. The procedure for making the DEM (Kirk et al., 2008) consists of resampling the images to remove optical distor- tions and then assembling the segments obtained by the 10 sepa- rate red-band detectors of the HiRISE camera in the USGS software system ISIS 3. We then use the commercial stereoanalysis 270 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 Fig. 1. Examples of light and dark striping in the NPLD at different scales, (a) Viking image 062B32. Image width = 68 km. Small, black box shows location of MOC image in lb. The residual cap is also called the Planum Boreum 4 unit in Tanaka et al. (2008); the NPLD on the south-facing trough walls in this location are part of the Planum Boreum 1 unit, while some of the north-facing trough walls expose the Planum Boreum 3 unit, (b) Portion of MOC image MOO/02072, showing NPLD layers in the Planum Boreum 1 unit. Image width = 0.83 km. Illumination is from the upper right. software, SOCET SET (?BAE Systems), to: control the images to Mars Orbiter Laser Altimeter (MOLA) data, produce a preliminary DEM by automated image matching, and perform interactive qual- ity control and editing of the DEM to remove matching errors. We use a grid spacing of 1 m between height points (posts) for the out- put product, and topographic features larger than a few meters across are fully resolved. The vertical precision (i.e., RMS amplitude of random errors in the DEM) is ~0.3 m. Taking into account the partial correlation of DEM errors between adjacent posts (Kirk et al., 2003), the precision with which short-baseline slopes can be estimated from the DEM is 2-3?. 3. What is a layer? For simplicity, in this paper, the NPLD refer to the ice-rich, lay- ered materials of Planum Boreum lying above the dark basal unit (Byrne and Murray, 2002; Fishbaugh and Head, 2005); all NPLD exposures that we discuss here are comprised of Tanaka et al.'s (2008) Planum Boreum 1 unit; see Tanaka et al. (2008) for a more detailed breakdown of the NPLD and basal unit into more specific geologic units. The martian polar deposits constitute the largest surface water reservoirs on the planet and thus are intimately linked to climate and hydrological cycles. The polar layers are likely to have re- corded the recent climate history, but before one can begin to understand that record, layers must be delineated. There is no sin- gle type of layer in a terrestrial ice sheet (e.g., Reeh et al., 1991; Svensson et al., 2005; Fischer et al., 2007; Lambert et al., 2008). An- nual snow accumulation, melting and re-freezing events, deposi- tion of dust and volcanic ash, and other processes all result in layering, some of which can be detected visibly, some by chemical analysis, some by radar, and so on. Because not all layers are visi- ble, we must do away with the preconceived notion that the dark and light stripes seen in Viking and Mars Orbiter Camera (MOC) images are the only "layers" in the martian polar layered deposits. Scale is also important. Many layers revealed in Viking images are actually packages of MOC-scale layers. If one used a method similar to Laskar et al. (2002) and attempted to tie Viking-scale lay- ering (defined by brightness) directly to changing orbital parame- ters or insolation and then did the same for MOC-scale layering, the results could likely be quite different. Fig. 2. (a) HiRISE stereo DEM created from stereo pair PSP_001871_2670 (87.0895?N, 92.7255?E) and PSP_001738_2670 (87.0946?N, 92.8213?E) overlying shaded relief created from the same DEM. Images were taken during Northern summer, separated by 5.316? Ls. Labeled lines are locations of profiles in Table 1 and Figs. 15 and 16. (b) MOLA data of same area for resolution comparison, (c) Context map is portion of MOLA 512 pix/deg shaded relief, (d) Shaded relief map created from the HiRISE DEM (using ESRI's ArcGIS 9.2), artificially illuminated with the same conditions as image PSP_001871_2670 and useful for assessing surface roughness. Surprisingly, HiRISE has not revealed even thinner layering than has MOC, but rather has highlighted much more detail of layers of the same scale as visible (some just barely) in MOC; the thinner layering that probably exists is obscured by younger, surficial deposits of dust and ice/frost and by slumping of this material and possibly of material inherent to the layers themselves (Herkenhoff et al., 2007). Here we refer to these deposits as the younger mantling deposits (YMD), some of which may be com- posed of dark veneers described by Rodriguez et al. (2007). The presence of this material is subtle and is mostly inferred from the gradational nature of layer boundaries and the absence of layers thinner than those visible in MOC images. However, as K.E. Fishbaugh et at./Icarus 205 (2010) 269-282 271 Fig. 3. (a) Example of a marker bed whose appearance changes within and outside of a dark streak of YMD material. The arrow separates the area within the dark streak (on the left) from the area without. PSP_001488_2665, image is 100 m across, illumination is from the lower right, (b) Eolian ripples within the YMD in image PSP_001871_2670, oriented perpendicular to bedding strike. Illumination is from the lower right. Image width is 150 m. 0.4 Brightness Fig. 4. (a) Orthorectified HiRISE image PSP_001738_2670 in uncorrected I/F units. Taken during Northern summer (Ls = 153?, illumination is from the lower left). Numbered locations are described in the text, (b) An aspect map (azimuth direction that each pixel is facing) made from the high-resolution stereo topography shown in Fig. 2. The DEM was binned to reduce noise in this product. Direction is measured in degrees east of north (180 is south-facing, 90 is east-facing, etc), (c) A slope map constructed from the same dataset. (d) Cosine of the incidence angle based on the illumination conditions in (panel a) and the slope (panel c) and aspect (panel b) of each point, (e) Lambert albedo derived by dividing (a) by (d). (f) Histograms of DN values for panels a,d and e. 272 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 ?! ? ??'? '?'?' 1 o ? Fig. 5. Qualitative comparison between an image and a portion of the shaded relief map (shown in full in Fig. 2d) shows that generally smooth and flat areas retain/collect more frost/ice than rougher areas, (a) Close up in shaded relief map, with illumination from the right, (b) Portion of HiRISE image PSP_001871_2670. Upslope direction is to the left. Illumination is from the right. Image width = 2 km. discussed further below, the YMD is sometimes directly recogniz- able as the material comprising some dark streaks (Fig. 3a), as sur- ficial coatings of frost/ice (Fig. 4a and Fig. 5b), and as ripples on layer surfaces (Fig. 3b). Given the different scales of layering (packages of thin layers, massive beds, etc.) for our purposes, we can only define a "layer" as a stratum evident in images and topography that cannot be di- vided into thinner strata at the best data resolution available. It thus has a generic definition, not attached to any particular charac- teristics and not implying any particular depositional/erosional environment. Each type of layer or layer package will likely have its own history of deposition, ablation, and perhaps alteration. One should keep in mind that given the resolutions of the available datasets, the presence of the YMD, and the fact that the water con- tent of the atmosphere is measured in precipitable microns, we probably cannot see annual layers (and possibly not even decadal layers). In this study, we take a significant step toward unraveling the climate record within the NPLD by delineating visible layers, as de- fined above, by their morphology and topographic expression. For these visible layer delineations to be useful for relating them to contemporaneous climate, they must represent isochrones, whether formed by deposition, as a lag after erosion, or during widespread alteration of pre-existing material. As discussed below, several of the layers have been identified and correlated across wide areas of the NPLD by Fishbaugh and Hvidberg (2006) and are therefore certainly isochrones. Determining the processes that created each layer is beyond the scope of this study but is a neces- sary next step. 4. High resolution layer morphology Many previous authors have used apparent brightness com- bined with layer morphology to identify particular layers in MOC images and to correlate them from place to place (Howard et al., 1982; Fenton and Herkenhoff, 2000; Malin and Edgett, 2001; Kolb and Tanaka, 2001; Milkovich and Head, 2005; Fishbaugh and Hvid- berg, 2006; Milkovich et al., 2008). Fishbaugh and Hvidberg (2006) have shown that the morphology of a particular layer in a MOC im- age can vary from one location to another. Even the well-known marker bed (Malin and Edgett, 2001), so called because of its char- acteristic knobby morphology in three different MOC images, changes appearance dramatically across the NPLD in MOC (Fishb- augh and Hvidberg, 2006) and HiRISE images. Because of the vari- ation in morphology, even of a particular layer, Fishbaugh and Hvidberg (2006) concluded that apparent brightness, morphology, and stratigraphic position are needed to identify specific layers and to correlate them across the NPLD. What these various authors have identified as "layers" in MOC images likely consist of a mix of many kinds of layers that can be defined in many ways: single, massive depositional beds, packages of thinner morphological lay- ers, layers with and without chemical differences from the adja- cent layers, morphologic layers and packages of layers that correlate with radar layers, etc. Using HiRISE images, we have been able to discern more details of the MOC-scale, visible layering. 4.1. Controls on brightness Analysis of HiRISE images has confirmed that apparent layer brightness is not necessarily indicative of the bulk composition of the layer, partially due to the presence of the YMD (Herkenhoff et al., 2007). Here, we use the HiRISE DEM and one of the stereo images from which it was created to further investigate potential controls on apparent layer brightness. Throughout the following discussion, we assume that without any surface covering of YMD, the inherent, aspect-independent albedo of any given layer would remain relatively constant along the entire layer exposure. Focusing on the layers exposed in the image in Fig. 4a, apparent layer brightness generally decreases downslope. Apparent bright- ness also decreases along layer strike, toward the bend in the trough at point 2 (Fig. 4a). Thus, any single layer does not have the same apparent brightness along its length in one image, and trough wall shape appears to exert some control on brightness. Part of the control of trough wall aspect on apparent layer brightness could be due to illumination conditions. To examine the influence of illumination conditions on apparent brightness, it is necessary to convert the radiometrically corrected image from l/F to aspect-independent albedo (i.e., to remove the effects of topographic slope and aspect on brightness). For HiRISE imagery, I/F is defined as the brightness observed (by HiRISE) relative to K.E. Fishbaugh et at./Icarus 205 (2010) 269-282 273 the brightness of a Lambert surface at the same distance from the Sun illuminated at 0? incidence, where / is the observed intensity and KF is the solar irradiance at Mars' distance from the Sun fil- tered through a HiRISE spectral bandpass (e.g., the red color band) (Delamere et al., 2009). There are significant complications in this conversion, which include estimating an unknown surface phase function, the amount of light scattered from the atmosphere to the surface, and the amount of light scattered from adjacent topog- raphy to a point of the surface. Achieving this level of accuracy is beyond the scope of our current study. Instead we have chosen to simplify this problem by assuming a Lambertian surface phase function and no incident radiation scattered from the atmosphere or adjacent topography. In this simple case, the (Lambert) albedo of a surface element is given by dividing the //F value by the cosine of the incidence angle. We calculate the incidence angle for each pixel by utilizing the DEM, hence aspect and slope (Fig. 4b and c) are ta- ken into account. The stereo processing procedure also produces an orthorectified image (i.e., Fig. 4a, a simulated nadir view, which is exactly coregistered with the DEM), and we take our IjF values from this orthorectified dataset. The DEM has a spatial resolution of one meter and a vertical accuracy of ~0.3 m. Calculation of slopes tends to amplify this noise as these uncertainties add in quadrature when figuring an elevation difference. This noise is considerably exacerbated when calculating a map of cosine of the incidence angle from the slope map. To mitigate the effects of this topographic noise, we have binned the DEM to 6 m/pixel be- fore calculation of the map of cosine of the incidence shown in Fig. 4d. Reducing 36 DEM pixels to one and extending the baseline over which slopes are calculated by a factor of six reduced the gra- dient errors by a factor of 36 and yielded considerably less noisy maps of incidence angle. Before dividing the 1/F values from the othrorectified image by the cosine of the incidence angle we re- stored the pixel size of the latter dataset to 1 m/pixel. In short, the Lambert albedo map shown in Fig. 4e is a combination of 1 m/pix brightness data and 6 m/pix topographic data. A comparison of the brightness data (Fig. 4a) and derived albedo data (Fig. 4e) shows that some of the large-scale brightness varia- tions were indeed due to the specifics of the trough wall's shape and orientation. For example, the general decrease of layer albedo with depth is much less prominent, and the dark region labeled 3 (Fig. 4a) is much brighter in Lambert albedo, both due to the re- moval of illumination effects. However, the albedo variation between frost-covered (e.g., point 1, Fig. 4a) and unfrosted portions (point 2, Fig. 4a) of the scarp wall persists, though to a lesser degree. Since the derived al- bedo does not remain constant along single layers, such albedo variations are most likely a result of the distribution of surface frost/ice and dust, rather than being due to variations in inherent layer albedo. This distribution might also be indirectly influenced by trough wall shape through control on wind patterns. Not all of the variance in image brightness can be explained by surface frost distribution and illumination conditions. Quantita- tively, the Lambert surface explains only 13.5% of the variation in image brightness. Fig. 4f shows that the brightness of this scene is bimodal. If we consider only the defrosted terrain (I/F less than the mean I/F of the scene) then the Lambert surface explains 52.5% of the brightness variance. The remaining unexplained vari- ance in the defrosted areas could be due to variation in inherent layer albedo, but is more likely due to the poor choice of surface phase function and variations in roughness, as described below. Fig. 4a shows the upper section of the trough wall varies in sur- face brightness on a small scale. In Fig. 5, a close up of the shaded relief map and the image, brightness appears related to the general roughness of the layers, with smoother surfaces collecting/retain- ing more frost/ice than rougher surfaces; rougher surfaces may also appear darker due to shadows cast by that roughness. Fig. 6 3 350 300 250 o Qi a. V) < *" - ."^"'Vy"^W T== ^"^f^^^".?^^^. I' SE + S n SW ? W ? NW ? N -NEiE bso 40 , ? j* v**? * ??9-* 110 250 ON Value Fig. 6. Illustration of the lack of correlation between DN value and meter/pixel physical parameters (slope and aspect). DN values from HiRISE image PSP_001871_2670. (a) DN value vs. aspect for a profile taken down the trough wall, across many layers. Aspect calculated from the DEM using ESRI's ArcGIS 9.2; units are degrees from longitude of center of image, (b) DN value vs. pixel-to-pixel slope (roughness) for same profile as above, calculated from the DEM using ArcGIS 9.2. (c) DN value vs. slope and aspect along one layer (profile width is one pixel), rather than down the trough wall. explores the potential relationships between brightness and layer properties derived from the DEM on an even smaller scale, pixel- to-pixel (1 m). Despite noisiness in the data as described above, it is clear that brightness is poorly correlated with the pixel-scale physical characteristics of the layers. A range of DN values is asso- ciated with each aspect direction (Fig. 6a), and Fig. 6b shows that slopes tend to cluster around 11? correlation with DN value. In Fig. 6c, we have plotted the DN value, slope, and aspect along the entire length of one layer (rather than cutting across many layers, down the trough wall slope). Again, while there is a noticeable 274 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 change in DN value along the layer (which is roughly correlated with the overall, curving wall shape), pixel-to-pixel scale slope and aspect appear uncorrelated to brightness. While, as discussed above, trough shape and general roughness have considerable influence on the distribution of the YMD and hence on surface albedo, the inherent properties of the layers must also have an influence, since correlations using MOC images (Fishbaugh and Hvidberg, 2006; Milkovich and Head, 2005; Milko- vich et al., 2008) have shown some consistency in layer brightness from place-to-place. But it is difficult to determine whether this brightness is inherent to the layer or if it only indirectly reflects layer properties by their control on distribution of the YMD and surface frost. Additionally, thus far only one NPLD HiRISE DEM has been created, so analysis of new OEMs may lead to slightly modified results. 4.2. Layer contacts As is evident from Fig. 7, delineating layers is difficult. For the most part, layer boundaries are not sharp at this scale but rather show a relatively gradual transition or even some overlap (the lat- ter likely being due to the YMD). Since there is likely to be a time lag in the response of the NPLD to changing climate, especially if the transition from one layer to another results from a relatively small climate shift, the contact between layers might not be ex- pected to be sharp when viewed at high resolution. However, if erosion produces a dusty lag centimeters thick, the transition from the underlying ice-rich layers to the lag would be sharp. It is imperative to note that the surface expression, when viewed from above, of a thin, near horizontal layer on a gently sloping surface will be greater than its true thickness. Additionally, and less obvi- ous, the surface of a layer may be gently sloping, while its edge is steep, creating a situation wherein the gently sloping surface ap- pears in nadir view as one layer, and the edge appears as a another layer below it, as illustrated in Fig. 8. In order to definitively estab- lish layer thicknesses and separations, a DEM is needed, as de- scribed in Section 5. 4.3. The original marker bed and other "marker beds" Malin and Edgett (2001) discovered a layer that was easily traceable in three different MOC images along one wall of a trough and named it the "marker bed". Milkovich and Head (2005) and Fishbaugh and Hvidberg (2006) later used this bed in their NPLD-wide layer correlations. As noted by Fishbaugh and Hvidberg (2006) the marker bed, identified by its stratigraphic position, changes appearance in MOC imagery from one location to the next, its only defining characteristics being a relative resistance to ero- sion and low apparent albedo. To identify the marker bed in HiRISE images, we use the MOC layer correlations provided in Fishbaugh and Hvidberg (2006) for stratigraphic context. Those authors iden- tify two separate stratigraphic sequences, the Upper and Lower Layer Sequences, the former of which contains the original marker bed. Both of these layer sequences lie within the Planum Boreum 1 unit of Tanaka et al. (2008). Fig. 9 shows examples of the marker bed from different locations throughout the NPLD in HiRISE images. If this layer formed in the same way at each location, the difference in appearance is likely primarily due to differences in lo- cal amount of YMD, intensity of erosion, slope, and lighting angle. Fig. 9 illustrates the difficulty in identifying particular layers in dif- ferent locations by their morphology alone. However, layers exposed within one image can easily be com- pared to one another. A close look at image PSP_001488_2665 (Fig. 10) for example, reveals that there are several layers similar in appearance to the original marker bed (the knobby layer in ^sloping layer surface Jayer edge Fig. 7. Example of PLD layers in HiRISE enhanced color. Grayscale inset shows close-up outlined by box and illustrates the lack of sharp transitions between layers. HiRISE image PSP_001534_2730 (87.0194?N, 195.467?E). Image width = 1.2 km. Inset width = 440 m. Illumination is from the upper right. Fig. 8. (Top) Illustration of a nadir view of layers exposed in a trough wall with example changes in albedo (shaded) and roughness (stippled), (middle) Illustration of a layer cross-section that interprets gently sloping, smooth layer surfaces as separate beds from the steeper, rougher layer edge faces. Differences in interpre- tation can cause confusion when delineating layer boundaries and counting layers, (bottom) Alternative interpretation in which gentle and steeply sloping sections are interpreted as the surface and edge faces of the same layer. K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 275 . #-\, Fig. 9. Examples of the differing appearance of the original marker bed (discovered in one trough by Malin and Edgett (2001)) from different locations around the NPLD. (a) PSP_001488_2665 (86.5101?N, 80.0374?E). Image width = 100 m. Upslope direction is toward the top. Illumination is from the bottom, (b) PSP_001530_2670 (87.1021?N, 358.639?E). Image width = 200 m. Upslope direction is toward the top. Illumination is from the lower right, (c) PSP_001760_2670 (87.011 ?N, 195.32?E). Image width = 200 m. Upslope direction is toward the upper left. Illumination is from the bottom right, (d) PSP_001554_2665 (86.2945?N, 88.1969?E). Image width = 250 m. Illumination is from the lower left. Upslope direction is toward the top. the center of Fig. lb). These layers are hummocky, exhibit linear erosional fluting on their upper edges, and, compared to the sur- rounding layers, are generally smoother, protrusive, and covered with less ice and frost (generally have a relatively surface low albe- do). Indeed, we have identified ten such layers in that image alone, leading to the revelation that there is more than one marker bed. The similarity in morphological appearance between these layers, a more obvious similarity than amongst other layers, also suggests that the marker beds may all have formed in the same way. To place these layers in stratigraphic context, in Fig. 10 we denote each marker bed and label those identified in MOC images by Fishbaugh and Hvidberg (2006). HiRISE has allowed what was pre- viously impossible?to identify these marker beds as having similar appearance, at scales much smaller than their thicknesses, and thus likely having the same origin, albeit separated in time. An- other set of marker beds (Fig. 11), perhaps of the same origin as the upper marker beds though not necessarily so, lie stratigraphi- cally below those in Fig. 10, within the Lower Layer Sequence. Again, we label these layers according to Fishbaugh and Hvidberg (2006). We have not found any evidence thus far of visible finer- scale layering within the marker beds, indicating that either they were deposited quickly, as massive beds, or that the YMD is shrouding the evidence of finer strata. It is not clear how much of their appearance is a result of the deposition of YMD on top of the actual layers. For example, in Fig. 3a, a marker bed exhibits a marked change in brightness at the boundary between a dark streak and the background layers and even appears to be smoothed by the streak (i.e., the streak has filled in hummocks in the layer exposure). This dark streak may be the result of wind action on the YMD and/or movement of dark veneer deposits (Rodriguez et al., 2007). The hummocky appearance could also result from slumping of the YMD or of the NPLD themselves, though we have found no evidence of such processes actively occurring. We have also observed thinner layer sets (layers ~1 m or less in thickness) between the marker beds in both the Upper and Lower Layer Sequences (Fig. 12), termed "laminated layers" by Milkovich and Head (2006). Due to erosion and YMD cover, these layer sets are difficult to correlate from one location to another. However, as more HiRISE images are acquired, it may become possible to at least detect the presence of thin layer sets in similar strati- graphic positions across the NPLD. 4.4. Erosional styles Erosion of the NPLD manifests itself at various topographic scales: scarps, troughs, and reentrants; erosion of layers into and out of the plane of a trough or scarp wall, forming apparent wavi- ness in the layer (though some such features may be caused by folding); and pits (perhaps formed by sublimation of ice), knobs, and elongated grooves at the small scale. As recognized in MOC images by Milkovich and Head (2006), the prevalence and particu- lar form of the pits, knobs, and roughness varies somewhat from layer to layer, likely due to differences in erosional resistance be- tween the layers (see inset in Fig. 7, for example). However, based on preliminary analysis of a few HiRISE images, the position and orientation of a particular outcrop exerts a larger control on ero- sion. Wind ripples have formed on layers outcropping in the walls of the one of the Olympia Cavi, whereas pitting is more prevalent 276 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 Fig. 10. Examples of multiple marker beds in one image in the Upper Layer Sequence (PSP_001738_2670). Labels in black on the left are the layers identified by Fishbaugh and Hvidberg (2006) in MOC images but not recognized, then, as marker beds. Labels in white on the right are layers newly identified in this study. All layers marked with white arrows are marker beds. Letters on far right denote examples shown in close-ups below. (B) is the original marker bed discovered by Malin and Edgett, 2001. Full image is 6 km across; close-ups are each 170 m across. Illumination is from the lower right. where those same layers wrap around to the west on a shallower slope (Fig. 13). Pervasive cracking in the steepest sections of the lower NPLD (Herkenhoff et al? 2007) in some Olympia Cavi scarps and in Chasma Boreale has overprinted erosional patterns in these layers. However, wind scouring is evident on the top surface of the PLD leading to the Chasma Boreale scarps (Fig. 14). From these observations, we suggest that wind erosion has a stronger influ- ence in areas likely to be experiencing greater katabatic wind velocities (i.e. steep scarps) (e.g., Howard, 2000) and that polygonal cracking such as that exhibited by the lower NPLD is slope-con- trolled. Erosion of the YMD may also account for much of the small-scale, erosional morphology. K.E. Fishbaugh et at./Icarus 205 (2010) 269-282 277 Fig. 11. Examples of multiple marker beds in one image in the Lower Layer Sequence. These marker beds are not necessarily of the same origin as those in the Upper Layer Sequence (Fig. 10) but are of the same origin as each other. Labels on the left are layers identified by Fishbaugh and Hvidberg (2006) in MOC images. Letters on right denote examples shown in close-ups below. Enhanced color, PSP_001636_2670 (84.137?N, 259.701 ?E, Ls = 143.918). Image is 1.2 km across. Grayscale close-ups are 100 m across. Upslope direction is toward the upper right. Illumination direction is from the upper right. 5. Topographic expression The creation of DEMs from HiRISE stereo pairs has significantly improved our ability to determine layer thickness and more accu- rate layer elevations than was previously possible using MOLA data, enhancing our ability to identify particular layers in various outcrops. However, the topographic expression of morphologically defined layers is highly variable. Fig. 15 shows several profiles ta- ken across a small portion of the DEM. It is immediately clear that the topographic expression of most of the morphologically visible layers is on the scale of the (erosional) roughness (~^1 m). In other words, it is difficult to pick-out layers in the topographic profiles. To reduce some of this erosional noise, we have averaged five sets of two profiles each. To measure layer thicknesses, we look for breaks in slope that appear in nearly the same location in each averaged profile, checking these locations in the HiRISE image and in the shaded relief created from the DEM to make sure that they lie within the morphologic expression of particular layers. We have measured maximum thicknesses, taking the topmost ele- vation of each layer as its top boundary; in some cases, these mea- surements may thus be including YMD deposited on top of the actual layer. Analysis of these profiles shows that major breaks in slope that exist in nearly the same location in every profile corre- spond to particular morphologic/stratigraphic layers previously identified in MOC images (Fishbaugh and Hvidberg, 2006), most of which correspond to what we can now identify as marker beds, as discussed above. Our method of measuring layer thicknesses is illustrated in Fig. 16. Note that in Fig. 16, the layer exhibits a peaked edge, and we use these peak edges, wherever existent in any layers, to estimate maximum layer thickness. The origin of these peaks is unclear, but they may form by erosion of the top sur- face of the layer, leaving a raised edge, or possibly by buildup of dust at rough layer edges. It is difficult to envision a primary, depo- sitional origin for this morphology. We summarize the results of these thickness measurements in Table 1, reporting results to the nearest meter. Layer thickness can naturally vary along one trough wall exposure and may ac- count for some of the difference in the measured thickness of one layer in different profiles. Errors in the determination of the elevation of layer boundaries are human errors and thus non-deterministic (DEM vertical precision errors (?0.3 m) are less than the layer thickness values). To estimate this error, layer boundaries were denoted independently by another coauthor on the shaded relief map alone, as shown in Fig. 17. The layer ele- vations measured in this way agree to within ~1 m of the eleva- tions shown in Table 1, hence the estimated 1 m error in the table. Notice in Fig. 16 that the brightness of the original marker bed changes halfway down its steep margin slope, further indicating that brightness and morphology alone are not enough to delineate and define layers. Indeed, when we have attempted to measure layer thicknesses simply by marking the locations of their albe- do/morphologic boundaries as apparent in the image, ignoring topography, the results differ for some of the layers. Layer thick- nesses tend to be slightly overestimated by that method because of the morphologically diffuse layer boundaries and covering of YMD, underscoring the value of high-resolution topographic data in stratigraphic analysis of the PLD. The next step of analysis is to create more HiRISE stereo DEMs to determine whether this pattern of layer thicknesses and separa- tion distances repeats across the NPLD and whether patterns also exist deeper in the stratigraphy. One can also make comparisons of the elevations of these prominent morphologic layers with the 278 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 ?crilrf ' Fig. 12. Examples of relatively thin layers, (a) Example in the Upper Layer Sequence. PSP_001630_2730 (87.0548?N, 96.0091?E, Ls = 143.683). Upslope is toward the upper left. Illumination is from the lower right. Image width is 100 m. (b) Example from the Lower Layer Sequence. PSP_001454_2670 (87.0395?N, 280.946?E, Is = 136.881). Illumination direction is from the lower right. Upslope direction is towards the upper left. Image width is 150 m. ????I m # \<* .vtf &*" ?Chasma Boreale wall 5 km | lvfl$LA 5 km! i~~n;? shaded relief SIR ""I , Fig. 13. Example of change of erosion style within the same layers exposed in outcrops with different aspects and slope PSP_001628_2650 (85.099?N, 236.753?E, Ls = 143.604). Illumination is from the lower left. Image width is 300 m. Black arrow shows downslope direction of steeper exposure (along an Olympia Cavi wall). White arrow shows downslope direction of more shallow exposure. Context maps show location of image. elevations of radar layers. If the morphologic layers correspond to radar layers, then we can also determine something about their dielectric properties, in turn helping us to determine their origin and place in the climate history. Fig. 14. Linear features on surface of NPLD, likely indicative of wind action. Tanaka et al. (2008) have interpreted similar such features as yardangs. Arrow shows direction with which linear features are aligned, roughly perpendicular to the steep Chasma Boreale wall slope. PSP_001334_2645 (84.4113?N, 343.516?E, Ls = 132.324). Image width is 800 m. Illumination direction is from the upper right. 6. Connections with climate The relationship between stratigraphy and climate may be quite complicated. Several factors can contribute to this complexity: (1) erosion of the NPLD and later deposition of YMD, potentially eras- ing some of surface expression of thinner layers and their bound- aries, resulting in apparently undifferentiated sections; (2) the fact that layer deposition most likely depends on the period AND amplitude of the polar insolation curve and other orbital parame- ters, possibly in a non-linear way and with lags in response time; (3) (related to 2) the fact that particular peak values are rarely re- peated in the insolation curve of the past 20 Myrs (Laskar et al., 2004); (4) episodic events, such as impacts, melting, periods of fas- ter ice flow, and volcanic eruptions; and (5) changing size and loca- tion of planetary water reservoirs contributing to NPLD build-up. Proper treatment of the complex problem of quantitatively tying K.E. Fishbaugh et at./Icarus 205 (2010) 269-282 279 -2615 Distance = 100 m Fig. 15. Illustration of the difficulty of separating topographic expression of layers from erosional "noise" on the scale of ~1 vertical meter. Five adjacent elevation (m) profiles down a short section of layers. Along-strike distance between top and bottom profile is ~845 m. Profiles have been shifted on the plot along the x-axis for easier viewing. No vertical exaggeration. See Fig. 2 for labeled profile locations, though note that this figure exhibits only a short section of the profiles. observed stratigraphy to climate is beyond the scope of this paper but is the subject of ongoing study (e.g., Fishbaugh et al., 2009). Rather, here we present a few preliminary, somewhat speculative scenarios that may explain some of the observed stratigraphy and pave the way for further discussion. One clue to a connection of NPLD stratigraphy with climate is the presence of the disconformity between the Upper and Lower layer sequences observed by Fishbaugh and Hvidberg (2006), which represents a major erosive episode and/or a change in depo- sition style and pattern. It is worth noting that the Upper and Low- er Layer Sequences together correspond roughly to the topmost radar-layered section observed by SHARAD (Phillips et al., 2008); the NPLD below that section exhibit a relative dearth of radar lay- ering, perhaps indicating a change in deposition style. A potential scenario to explain these major unconformities is illustrated by Fig. 18a. Perhaps the marker beds are deposited during times of moderate obliquity (e.g., B and D), with their spacing controlled by the period of the obliquity cycle. In that case, perhaps the times of very small obliquity fluctuations (in which state Mars is at pres- ent, e.g., A and C) create the disconformities between layer se- quences. It should be noted that small-scale, angular unconformities do not appear to be widespread within north polar stratigraphy and are concentrated mostly along the Planum Bore- urn margins (Milkovich and Head, 2005; Tanaka, 2005; Fishbaugh and Hvidberg, 2006), and we do not consider their implications for climate history in this study. Alternatively, Cutts and Lewis (1982) model the deposition of icy, dusty layers under changing orbital parameters. In their Climate Modulated Deposition Model, the major layer constituent, ice, is controlled by climate changes, while the minor constituent, dust, has a constant accumulation rate. Layer boundaries mark a decrease in or lack of ice deposition. The model simply allows ice to ablate and accumulate based assumed accumulation rates and on ablation rates calculated by Toon et al. (1980), which are al- lowed to vary with the orbital parameters. The model predicts that layers will generally accumulate in a repeating pattern of a rela- tively thick layer followed by thinner layers. The thick layer forms when the amplitude of obliquity variations is small, presumably due to the low ablation rate (though accumulation is also relatively slow). The difference in thickness between the "thick" and "thin" layers is influenced primarily by the assumed accumulation rate; the higher the accumulation rate (i.e. the lower the climate-con- trolled threshold), the greater the dichotomy in layer thickness. Layers formed when strong precession control coincides with low amplitude obliquity variations tend to be thinner and more closely spaced. How do Cutts and Lewis (1982) model results compare with what we actually observe? We have indeed observed relatively thick beds (the marker beds) and intervening thinner beds. The marker beds could thus have formed during times of low ampli- tude obliquity variation. If we assume this to be true, then accord- ing to Fig. 18b, only eight of the 10 identified marker beds from the Upper Layer Sequence could have formed in the last 20 Myrs. The full thickness of the NPLD would, in that case, be much older than 20 Myrs since the layer sequence covered by the DEM in this study spans only the top few hundred meters. Levrard et al. (2007) use a global climate model to build icy lay- ers of the NPLD, taking into account theoretical shifts in water sup- ply reservoirs from the tropics to the high latitudes at different obliquity values and the creation of dust lags when ice sublimates. Their model produces distinct layers associated with the 51 kyr and 120 kyr insolation cycles. These layers range in thickness from 10 to 80 m, influenced by the 2.4 Myr and 1.3 Myr obliquity mod- ulation periods. These layer thicknesses are greater than the range in marker bed thicknesses reported in Table 1, and the number of layers produced in the model is fewer than the number of layers 150 m Fig. 16. Illustration of the method used to measure layer thicknesses. This example shows the original marker bed. On the left are five averaged topographic profiles over the marker bed with gray circles denoting the marker bed boundaries as identified in HiRISE image PSP_001871_2670 and in the shaded relief map and by searching for breaks in slope in the topographic profiles. The along-strike distance between the two farthest profiles is 845 m (profile locations are marked in Fig. 2). Profiles have been shifted on the plot along the x-axis for easier viewing. On the right is a perspective view of the marker bed as exposed in PSP_001871_2670, created using the HiRISE DEM with a vertical exaggeration of 5x. The circles denote the marker bed boundary as identified above. Note that the top boundary of the marker bed has been marked at its peak elevation to gain a maximum thickness estimate. 280 K.E. Fishbaugh et aL/Icarus 205 (2010) 269-282 Table 1 Elevations (m), maximum thicknesses (m), and separation distances (m) of marker bed layers measured in the HiRlSE DEM. Layers E, Ul, and MB were identified in Fishbaugh and Hvidberg (2006) (see Figs. 10 and 11). Profile positions denoted in Fig. 2. Layers E Ul MB MB-1 MB-2 MB-3 Profile A Top Bottom Thickness -2606 -2613 7 -2643 -2650 7 -2681 -2689 8 -2712 -2716 4 -2732 -2738 6 -2744 -2750 6 Profile E Top Bottom Thickness -2608 -2616 8 -2644 -2650 6 -2682 -2690 8 -2713 -2710 3 -2729 -2733 4 -2740 -2746 6 Profile B Top Bottom Thickness -2607 -2615 8 -2643 -2649 6 -2683 -2690 7 -2714 -2718 4 -2731 -2735 4 -2746 -2750 4 Profile D Top Bottom Thickness -2608 -2617 9 -2645 -2649 4 -2682 -2690 8 -2714 -2710 4 -2734 -2740 6 -2742 -2746 4 Profile C Top Bottom Thickness -2608 -2617 9 -2646 -2652 6 -2684 -2691 7 -2715 -2729 4 -2736 -2741 5 -2740 -2746 6 Error(?) Top Bottom Thickness 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Average thickness Error 8 V2 6 V2 8 V2 4 V2 5 V2 5 V2 Separation distance Error 29 32 24 16 5 observed. According to Levrard et al. (2007), the latter discrepancy may be explained by orbital parameters other than obliquity (the main control on polar insolation) having a marked influence on layer deposition (as noted by Cutts and Lewis (1982)). Both prob- lems may also be explained by the complicating factors listed at the beginning of this section. Additionally, it is possible that such climate modeling explains the creation not of the actual, physical, individual layers, but rather of distinct layer packages, not as easily recognized as are individual layers. Spectral analysis of changing layer brightness with depth shows dominate depth wavelengths of ~30 m (Milkovich and Head, 2005; Milkovich et al., 2008) and ~20-80 m (Perron and Huybers, 2008) that may correspond more directly to the layering predicted by Levrard et al. (2007). In any case, the current inability of climate models to accurately predict the observed sequence of layers makes it clear that NPLD stratigra- phy does not result from a simple relationship between insolation and accumulation/ablation rate alone. The erosional resistance of the marker beds might also provide a climate connect clue. Milkovich and Head (2005) have proposed that sublimation of the NPLD during a period of net ablation re- sulted in the accumulation of a dark dust lag. If the relatively dark, erosionally resistant marker beds correspond to such a lag and if the NPLD contain only a few percent dust on average (Phillips et al., 2008), then ~400 m of NPLD material must ablate to form the thinnest marker bed in Table 1. Such massive amounts of ero- sion would need to occur periodically to form all of the marker beds. Also, since the marker beds do not cross-cut the layers below, downwards ablation would need to be extremely uniform laterally across the NPLD, making this scenario unlikely. Alternatively, the marker beds may have been created when the dust content of the atmosphere was higher (during higher obliquities), resulting in a higher dust/ice ratio in the PLD layers (but not much higher since the dust percentage overall is quite low) and a higher resis- tance to sublimation. Even more complicated scenarios could in- volve such factors as increased temperatures (again, during higher obliquities) resulting in ice metamorphism, creating larger ice crystals, which would decrease the albedo and perhaps anneal the grains (under partial melting conditions), strengthening the layer. Note that preferential accumulation of YMD dust and shad- ows created by surface roughness could also account for the low al- bedo of the marker beds. All of these scenarios involve formation of marker beds during high obliquity excursions, rather than during periods of low amplitude variation, as discussed above and illus- trated in Fig. 18b. So how and when does a thick, erosionally resistant layer form? The direct and definitive connection between deposition of particular layer types with particular thicknesses remains enig- matic and will be the subject of future study (e.g., Fishbaugh et al. (2009)). Speculatively, the separation distance between layers may be controlled by a complex interplay of several orbital parameter variations, but the thickness and physical properties of the layers may result from the particular amplitudes of these variations at the time when the layer is deposited. Note that, according to Laskar's calculations (Laskar et al., 2004), the ampli- tude of the insolation curve is not always the same at the two ends of one cycle. Recently, Lewis et al. (2008) observed repeating layer sequences in Arabia Terra and measured their thicknesses using a HiRlSE ste- reo DTM. Layer thicknesses there are typically ~10m, similar in scale to the polar marker beds, though in some places there are bundles of thick and thin layers with a 10:1 thickness ratio. Future comparisons between the two types of layering using more DEMs could prove useful if the deposition/erosion rates of both are con- trolled in large part by orbital parameters. 7. Summary 1. Layers are delineated by the method used to detect them. In this study, the PLD layers we delineate are strata that are evident in images and topography and cannot be subdivided into thinner strata at the best data resolution available. This definition does not exclude the possibility that thinner, undetectable visible strata exist in the NPLD or that other types of layering exist. 2. We find that not all of the variation in layer brightness can be explained by the inherent albedo of the layers. Indeed, most K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 a 50 281 Fig. 17. Illustration of selecting layer boundaries along a cross-section in a shaded relief map derived from the HiRISE DEM, using lateral continuity and roughness variations relative to adjacent layers. Shaded-relief derived from the DEM. of the brightness variation appears to be due to illumination angle (affected by trough shape and orientation) and the pres- ence of the YMD. The distribution of the YMD is likely in turn influenced by (1) the shape of the trough wall, due to its effects on trough-local wind patterns, (2) by the layer surface rough- ness, partially controlled by the inherent erosional resistance of the layers (note, however, that pixel-scale aspect and slope do not appear correlated with pixel brightness), and (3) possibly by such inherent layer properties as thermal inertia. Layer delineations cannot be based on morphology or brightness alone but require topographic data, indicating major control of surface YMD cover on layer appearance. C 40 ?5 30 a 20 o 10 Mil ,i 1 -4 -6 Time (Myr) -4 -6 Time (Myr) Fig. 18. Illustration of possible timing of marker bed formation in different scenarios. Obliquity values are taken from Laskar et al? 2002. Letters correspond to time periods discussed in the text, (a) One possible scenario in which marker beds form during periods of moderate obliquity (B,D): black arrow shows example of when a marker bed might form, and grey arrow shows an example of when the intervening layers might form. Major unconformities might be created during periods A and C, as, for example, between the Upper and Lower Layer Sequences (Fishbaugh and Hvidberg, 2006). (b) Another possible scenario in which marker beds form during periods of small obliquity variation (circles) and low orbital inclination values (see Laskar et al? 2004 for inclination calculations), as discussed by Cutts and Lewis (1982). 3. The Upper Layer Sequence of the NPLD contains several marker beds similar in appearance and likely in origin to the original marker bed discovered by Malin and Edgett (2001). Specula- tively, we suggest that these marker beds may correlate with the radar layering observed in the upper portion of the NPLD (Phillips et al., 2008; Putzig et al., 2008). From the DEM, we obtain marker bed thicknesses ranging from 4 to 8 ? \/2m and separation distances ranging from 5 to 32 ? A/2 m. Marker beds not necessarily of the same origin as those above them have also been discovered in the Lower Layer Sequence. 4. Layer contacts appear gradational in images and subtle in topo- graphic expression. This blurring of layer expressions is at least partly due to the presence of the YMD. 5. Preliminary observations indicate that small-scale erosion mor- phology (a mix of the morphology of the layers themselves and of the YMD) is qualitatively, largely controlled by the location within the PLD and orientation of the exposing trough wall or scarp. Scarps and walls that are likely exposed to greater kata- batic wind velocities exhibit more wind erosion than others scarps and walls that exhibit more (sublimation?) pitting. Com- prehensive mapping of erosional styles constitutes a detailed study beyond the scope of this paper. 6. A definitive connection between orbital parameter values and the deposition of particular layer types with particular thick- nesses cannot be established in a simple way. In other words, we find it unlikely that the period of any one orbital parameter, including obliquity, could be directly tied with DN value, layer thickness, or layer separation distance. It is likely that the latter two values are influenced by a climatic signal whose period depends on the combination of the periods of several orbital parameters and whose amplitude depends on the particular amplitudes of these orbital parameters at the time of layer deposition. 282 K.E. Fishbaugh et al./Icarus 205 (2010) 269-282 Acknowledgments The authors profusely thank those at the U.S. Geological Survey in Flagstaff who were involved in the creation of the stereo DEM: Elpitha Howington-Kraus, Donna M. Galuszka, Bonnie L Redding, and Mark R. Rosiek. Thank you also to the HiRISE team for so many helpful discussions. Two anonymous reviews greatly strengthened the paper. This work was funded by a NASA MRO Participating Scientist grant to KEF. References Byrne. S., Murray. B., 2002. North polar stratigraphy and the paleo-erg of Mars. J. Geophys. Res. 107 (E6). doi:10.1029/2001JE001615. Cutts, J., Blasius, IC, Briggs, G., Carr, M.. Greeley, R., Masursky, H.. 1976. North polar region of Mars: Imaging results from Viking 2. Science 194, 1329-1337. Cutts, J.. Lewis, B.. 1982. Models of climate cycles recorded in martian polar layered deposits. Icarus 50, 216-244. Delamere, WAand 15 colleagues, 2009. Color imaging of Mars by the High Resolution Imaging Science Experiment (HiRISE). Icarus 205, 38-52. Fenton, L, Herkenhoff, IC, 2000. Topography and stratigraphy of the martian north polar layered deposits using photoclinometry, stereogrammetry, and MOLA altimetry. Icarus 147. 433-443. Fischer. H., Siggaard-Andersen, M., Ruth, LI., Roethlisberger, R.. Wolff, E., 2007. Glacial/Interglacial changes in mineral dust and sea-salt records in polar ice cores: Sources, transport, and deposition. Rev. Geophys. 45 (RG1002). doi:10.1029/2005RG000192. Fishbaugh,!(., Head. J., 2005. Origin and characteristics of the Mars north polar basal unit and implications for polar geologic history. Icarus 174, 444-474. Fishbaugh,!(., Hvidberg, C, 2006. Martian north polar layered deposits stratigraphy: Implications for accumulation rates and flow. J. Geophys. Res. Ill (E06012). doi:l 0.1029/2005JE002571. Fishbaugh,!(., Hvidberg, C, Byrne, S., Herkenhoff, K., Fortezzo, C, Kirk, R., Winstrup, M., 2009. The stratigraphic record in the martian north polar layered deposits as measured by high resolution stereo topography. Lunar Planet. Sci. 40. Abstract #1998. Herkenhoff, K., Murray, B., 1990. Color and albedo of the south polar layered deposits on Mars. J. Geophys. Res. 95, 14511-14529. Herkenhoff, IC, Byrne, S., Russell, P., Fishbaugh, IC, McEwen, A., 2007. Meter-scale morphology of the north polar region of Mars. Science 317, 1711-1715. Howard, A.. 2000. The role of eolian processes in forming surface features of the martian polar layered deposits. Icarus 144, 267-288. Howard, A., Cutts, J., Blasius, IC, 1982. Stratigraphic relationships within the martian polar cap deposits. Icarus 50, 161-215. Kirk, R., Howington-Kraus, E., Redding, B., Galuszka, D., Hare, T., Archinal, B., Soderblom, L, Barrett, J., 2003. High-resolution topomapping of candidate MER landing sites with Mars Orbiter Camera Narrow-Angle images. J. Geophys. Res. 108 (E12). doi:10.1029/2003JE002131. Kirk, R., and 19 colleagues, 2008. Ultrahigh resolution topographic mapping of Mars with MRO HiRISE stereo images: Meter-scale slopes of candidate Phoenix landing sites. J. Geophys. Res. 113 (E00A24), doi: 10.1029/2007JE003000. Kolb, E., Tanaka, IC, 2001. Geologic history of the polar regions of Mars based on Mars Global Surveyor data. II. Amazonian period. Icarus 154, 22-39. Lambert, F., Delmonte, B., Petit, J., Bigler, M., Kaufmann, P.. Hutterli, M., Stocker, T., Ruth, U., Steffensen, J., Maggi, V., 2008. Dust-climate coupling over the past 800,000 years from the EPICA Dome C ice core. Nature 452. doi:10.1038/ nature06763. Laskar.J.. Levrard, B., Mustard. F., 2002. Orbital forcing of the martian polar layered deposits. Nature 419, 375-377. Laskar.J., Correia, A., Gastineau, M.,Joutel. F., Levard, B.. Robutel, P.. 2004. Long term evolution and chaotic diffusion of the insolation quantities on Mars. Icarus 170, 343-364. Levrard, B., Forget, F., Montmessin, F.. Laskar. J., 2007. Recent formation and evolution of northern martian polar layered deposits as inferred from a Global Climate Model. J. Geophys. Res. 112 (E06012). doi: 10.1029/2006/JE002772. Lewis, IC, Aharonson, O., Grotzinger, J., Kirk, R., McEwen, A.. Suer, T.-A., 2008. Quasi- periodic bedding in the sedimentary rock record of Mars. Science 322,1532-1535. Malin, M., Edgett, IC, 2001. Mars Global Surveyor Mars Orbiter Camera: Interplanetary cruise though primary mission. J. Geophys. Res. 106, 23429-23570. McEwen, A., and 14 colleagues, 2007. Mars Reconnaissance Orbiter's High Resolution Imaging Science Experiment (HiRISE). J. Geophys. Res. 112 (E05S02). doi:10.1029/2005JE002605. Milkovich, S.. Head, J., 2005. North polar cap of Mars: Polar layered deposit characterization and identification of a fundamental climate signal. J. Geophys. Res. 110 (EOS). doi:10.1029/2004JE002349. Milkovich, S., Head, J., 2006. Surface textures of Mars' north polar layered deposits: A framework for interpretation and future exploration. Mars 2. 21-45. S. Milkovich, J., Head, G., Neukum and HRSC Co-Investigator Team, 2008. Stratigraphic analysis of the northern polar layered deposits of Mars: Implications for recent climate history. Planet. Space Sci. 56, 266-288. Murray, B., Ward, W., Yeung, S., 1973. Periodic insolation variations on Mars. Science 180, 638-640. Perron, J., Huybers, P., 2008. Is there an orbital signal in the polar layered deposits on Mars? Paper presented at the Lunar Planet. Sci. 39 (Abstract #1391). Phillips, R., and 26 colleagues, 2008. Mars North Polar Deposits: Stratigraphy, age, and geodynamical response. Science 320, 1182-1185. Putzig, N., and 18 colleagues. 2008. Stratigraphic Mapping of the North polar layered deposits on Mars from Radar Soundings. Lunar Planet. Sci. 39 (Abstract #3295). Reeh, N? Oerter, H., Letreguilly, A., Miller, H., Hubberten, H.-W., 1991. A new, detailed ice-age oxygen-18 record from the ice-sheet margin in central West Greenland. Global Planet. Change 4 (4), 373-383. Rodriguez, J.. Tanaka, IC, Langevin, Y., Bourke, M., Kargel, J.. Christensen, P., Sasaki. S., 2007. Recent eolian erosion and deposition in the north polar plateau of Mars. Mars 3, 29-41. Svensson, A., Nielsen, S.. Kipfstuhl, S., Johnsen, S.. Steffensen. J.. Bigler, M., Ruth, U., Rothlisberger. R., 2005. Visual Stratigraphy of the North Greenland Ice Core Project (NorthGRIP) ice core during the last glacial period. J. Geophys. Res. 110 (D02108). doi:10.1029/2004JD005134. Tanaka, IC, 2005. Geology and insolation-driven climatic history of Amazonian north polar materials on Mars. Nature 437, 991-994. Tanaka, IC, Rodriguez, J.. Skinner, J., Bourke, M., Fortezzo, C. Herkenhoff. IC. Kolb, E., Okubo, C, 2008. North polar region of Mars: Advances in stratigraphy, structure, and erosional modification. Icarus 196 (2), 305-317. Toon, O.. Pollack, J.. Ward, W., Burns, J., Bilski, IC, 1980. The astronomical theory of climatic change on Mars. Icarus 44. 552-607.