1 Introduction

Carbon dioxide ice is observed on Mars in various locations and timescales, from ephemeral clouds, to diurnal frost and seasonal caps (Calvin & Titus, 2008), to the south polar residual cap (SPRC, Figure 1) (Byrne, 2009; Kieffer, 1979) and massive CO₂ ice deposits (MCID, Phillips et al., 2011). Most of these can be seen from Earth when geometry and seasons permit, but the thin, high albedo SPRC that persists over decadal time scales was of primary interest for early Mars observations because it may have implied that the entire south polar layered deposits (SPLD, Figure 1) were composed primarily of carbon dioxide ice. This led to laboratory investigations considering the rheological properties of CO₂ ice and its effect on the SPLD.

These investigations demonstrated that CO₂ ice may flow up to 100x faster than water ice at Martian temperatures (Clark & Mullin, 1976; Durham et al., 1999). Based on those strength measurements, Nye et al. (2000) ruled out the possibility of an SPLD composed entirely of CO₂ ice because it would have insufficient strength to maintain its observed shape over the long time periods implied by the SPLD's surface cratering age (Herkenhoff & Plaut, 2000). That study concluded that the CO₂ ice seen in the SPRC must be thin and cover a more rigid water ice cap, a conclusion later supported by higher resolution observations (Byrne & Ingersoll, 2003) and gravity observations that are sensitive to the entire SPLD thickness (Wieczorek, 2008; Zuber et al., 2007). The possibility of CO₂ glaciers was reconsidered when evidence for moraines near the north polar layered deposits was discovered (Kreslavsky & Head, 2011), but no CO₂ glaciers are present in the north today.

That was the state of knowledge until 2011, when observations by the Shallow Radar (SHARAD) instrument (Seu et al., 2007) on Mars Reconnaissance Orbiter (MRO) detected MCID that reside in deep, elongated, curvilinear basins and are buried beneath a thin layer of dusty water ice and the SPRC (Phillips et al., 2011, Figure 2). Using geophysical arguments and layer geometry, the real part of the dielectric constant of the material was determined to be consistent with that of CO₂ ice, up to 1,000 m thick locally but averaging ∼400 m thick (Bierson et al., 2016; Phillips et al., 2011; Putzig et al., 2018), two orders of magnitude thicker than the SPRC (Byrne, 2009; Byrne & Ingersoll, 2003). The thickness of these deposits varies in relationship to the depth of the basins in which they reside, and several isolated deposits are detached from the main, near-contiguous body of ice (Figure 2). The ice volume, as determined by measurements using a recently generated three dimensional (3D) SHARAD data set (Foss et al., 2017) and extrapolation toward the pole, exceeds 16,500 km³, equivalent to ∼2.7 × 10¹⁵ kg (Alwarda & Smith, 2021; Putzig et al., 2018). Using the 1,600 kg/m³ density of CO₂ ice (Piqueux et al., 2003), this mass is approximately equal to the mass of the current atmosphere (Guo et al., 2010).

Based on these detections, modeling efforts have attempted to explain the depositional and climatic history of the MCID (Buhler et al., 2019; Manning et al., 2019). They predicted similar timeframes for the initial deposition, which would have began ∼510 kyr in the past, but different mechanisms of deposition and preservation. Buhler et al. (2019) and a follow up study (Buhler & Piqueux, 2021) were able to calculate deposition rates through time for these deposits and a total mass that closely approximates that observed by radar (Putzig et al., 2018). The modeling of Manning et al. (2019) predicted that the bounding layers that separate each CO₂ unit would act as seals to protect the CO₂ ice during warm periods, similar to suggestions by Bierson et al. (2016); however, this came with a stratigraphic section of CO₂ ice that was sandwiched between thin H₂O layers. Alwarda and Smith (2021) compared the two models to SHARAD observations of the stratigraphy and found evidence to support the model presented in Buhler et al. (2019) but not that of Manning et al. (2019). Neither one-dimensional model treated the three-dimensional distribution of CO₂ ice that is observed, leaving that question open for future work. Work from a decade earlier suggested that the distribution of the SPRC, namely the offset from the geographical pole, the average elevation much lower than the highest elevation on the cap, and the presence between 120°W and 45°E (Figures 1 and 2) are caused by large scale topographic forcing of atmospheric dynamics (Colaprete et al., 2005). Extrapolating backwards in time and assuming a similar mechanism for the MCID and the SPRC, this explanation could explain the general location of the MCID but not the three-dimensional distribution, including the thickness or disconnected nature.

Here, through observations (Sections 2 and 3) and modeling (Sections 4 and 5), we demonstrate for the first time that the distribution of the MCID is caused by glacial flow and that uniform atmospheric deposition is insufficient to explain the observations. In Section 6 we discuss the implications of this finding and answer some previously outstanding questions about the distribution of the ice.

2 Geomorphic Observations

We employ datasets from ongoing and past Mars orbiters to characterize the physical state of the CO₂ deposits. Topography is measured with the 512 pixels per degree south polar data set from the Mars Orbiter Laser Altimeter (MOLA, Smith et al., 2001) and supplemented with a south polar digital terrain model (DTM) from the High Resolution Stereo Camera (HRSC, Neukum & Jaumann, 2004; Putri et al., 2019). Surface features are identified and interpreted with data from the High-Resolution Imaging Science Experiment (HiRISE, McEwen et al., 2007) and the Context Imager (CTX, Malin et al., 2007) on MRO. Subsurface geometry is measured with three-dimensional data (Foss et al., 2017; Putzig et al., 2018) from SHARD (Seu et al., 2007).

Using MOLA and orbital imagery, we find that the MCID exhibit variable physical characteristics that can be classified as laterally bound, laterally unbound, or partially bound. The first landform is a ∼100 m thick deposit that resides on the highest elevation and latitude of the SPLD plateau. Here, the CO₂ ice resides on a level surface and is unbound at the margins (Red Star #1 in Figure 3). Topographic measurements only cover a portion of this deposit, but where it can be measured, the gross morphology resembles a dome. The second type of morphologic expression, laterally bound deposits, is found in two locations at the lowest elevations within closed basins. The two examples express low relief and are physically disconnected from nearby deposits (Red Stars #2 and #3 in Figure 3). Third, at intermediate latitudes and elevations, several other deposits reside in elongated, curvilinear basins. These deposits are bound on two sides by preexisting curvilinear basins that resemble spiral troughs when seen in 3D radar (Alwarda & Smith, 2021), and the surfaces dip in the same direction as the underlying topography (Figure 2; Red Stars #4 and #5 in Figure 3). Further, each unit has convex up profiles, at least for the majority of their length, that point in the dip direction. We also note a concave portion in the highest 30 km of Profile 4. In some cases, the surface slope steepens slightly at the lowermost section (Figure 3, Stars #4 and #5). We informally name two of these partially bound deposits Wolf Glacier and Brooks Glacier.

Other features are unique to the thick CO₂ units as well. Over much of the SPRC there are shallow elongated depressions (minor troughs; Smith et al., 2014) that orient perpendicular to the maximum topographic gradient and parallel to the large-scale, parabolic topographic contours as mapped by MOLA (Figure 3 and Figures S1–S3 in Supporting Information S1). They are best seen in late summer as dark curvilinear marks within the high albedo SPRC (Figures S1 and S2 in Supporting Information S1), have depths that range from 10 to 100 m, and are hypothesized to be sublimation pits (Phillips et al., 2011) or erosional megadunes (Smith et al., 2014). At the highest elevations, near to the lowest slopes, the minor troughs have very small topographic perturbations (<10 m) and have short wavelengths, whereas the deposits at intermediate elevations and greater surface slopes exhibit much larger topographic depressions up to 100 m deep and tens of km wide (Profile #4 in Figure 3). The lowest elevation deposits that reside in closed basins exhibit no minor troughs.

Additionally, smaller surface features include low-relief ridges that vary from 5 to 20 m in height (Figure 4). Near the margins of Wolf and Brooks Glaciers, we observe low ridges that align parallel to the large-scale topographic contours (Figure 4). Similar features are also observed on the margins of the two isolated, lower elevation deposits (red arrows in Figure S4 of Supporting Information S1). The locations are indicated by yellow stars in Figure 3.

3 Geologic Interpretation

We interpret the physical attributes observed as indicative of current or relict flow of CO₂ glaciers. For the gross morphology, the landforms we observe are end members of viscous flow with different boundary conditions. The filled basins are the result of glaciers flowing downhill into basins where they are confined and will eventually become fully level. The dome-shaped deposit at the highest plateau of the SPLD resembles unconfined flow of CO₂ over a flat base (e.g., Nye et al., 2000). The elongated deposits at intermediate elevations resemble terrestrial outlet or valley glaciers (Figure S5 in Supporting Information S1) that have parabolic topographic contours pointing toward the glacier toe (Figure 3a and Figures S1–S3 in Supporting Information S1).

For the elongated, valley-filling deposits, the parabolic contours and the convex-up shape (Figure 3) are a result of the fluid responding to gravitational stress. The convex-up shape is indicative of an ablation zone in terrestrial outlet glaciers, whereas the concave-up region, near the head of Wolf Glacier, is more typical of an accumulation zone (Cuffey & Paterson, 2010, Figures 4.4 and 11.7). We don't believe that the glacier is accumulating from atmospheric deposition at present, so accumulation may arrive from further south at the ice sheet margin. Unfortunately, that is right at the latitude where MOLA and SHARAD lose coverage; however, topographic data from a HRSC DTM (Putri et al., 2019) suggests that the topography increases with latitude, potentially supporting this interpretation. Furthermore, for Wolf and Brooks Glaciers, we interpret the increased slopes at the terminus of the deposits as glacial toes that steepen because the base does not slide (Stars #4 and #5 in Figure 3).

We interpret the minor troughs (Figures S1 and S2 in Supporting Information S1) to be crevasses that form due to the strong topographic gradients that cause elevated stresses and fracturing. We consider flow direction (as demonstrated by the parabolic contours and our modeling results; see Section 5) and find that the crevasses resemble ogives at outlet glaciers in ablation zones (Figure S5 in Supporting Information S1); however we believe that the formation mechanism differs from that of ogives because ogives can only form near icefalls (Cuffey & Paterson, 2010), and that is not possible for the minor troughs. Fracturing occurs as a result of gravitational driving stresses overcoming the strength of a mechanically strong but thin (∼8 m thick) water-ice layer that overlies and protects the CO₂ deposits (BL₃ in Buhler et al. [2019] and Alwarda and Smith [2021]). Driving stresses are positively influenced by increased gravitational acceleration, surface slope, ice thickness, and ice density and temperature. Because the stresses point downhill, and the water-ice layer is prone to brittle failure at 150 K, the failures will be parallel to the topographic contours, as observed, matching crevasse formation on terrestrial valley glaciers (Figures S2 and S5 in Supporting Information S1). For more discussion about crevasse growth and alternative interpretations regarding the minor troughs, including winds and preferential deposition, see Supporting Information S1.

We interpret the low-relief ridges (Figure 4 and Figure S4 in Supporting Information S1) to be compressional features near the margin of the glacier. They resemble compression ridges found near McMurdo Station, Antarctica (Figure S6 in Supporting Information S1), where faster moving ice abuts slower moving ice. This type of buildup of stress usually happens at the lateral margin of a glacier or at the glacial toe in locations where sliding is not prevalent (Cuffey & Paterson, 2010, Figure 11.7).

In an attempt to detect flowing ice, we compared available HiRISE and CTX images of the CO₂ deposits for the longest baseline possible. Our search for surface translation over the relatively short baseline of Mars Reconnaissance Orbiter's time in orbit (∼15 years) has not detected any movement. This is unsurprising, as our modeling (See Section 5) suggests that maximum velocities today are <0.2 m yr⁻¹, for a maximum translation of <3 m over the lifetime of MRO. HiRISE has resolution capable of seeing this motion translation, but this is not a conservative estimate of flow, so it is likely slower, and we have not identified any noticeable translation yet.

4 Glacial Modeling Methods

To test the formation history of the MCID and SPRC, we perform three-dimensional glacial modeling using the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012). ISSM is a state-of-the-art finite element thermo-mechanical ice flow model, capable of simulating transient flow on a regional or continental scale. The ISSM software is open-source (http://issm.jpl.nasa.gov) and is designed as a flexible and scalable platform for the physical modeling of ice flow and the incorporation of data assimilation techniques via inverse control methods (Morlighem et al., 2009). ISSM has an adaptable, anisotropic mesh that can be highly refined in areas of complex ice flow or steep surface topography, and the software itself can handle simulations with millions of degrees of freedom, relying upon a C++ core for dynamic memory allocation and parallelization through multi-processor distributed computing. Ice flow can be modeled in 2 or 3 dimensions (2D/3D) with a number of momentum balance models including a Shallow-Shelf Approximation and the 3D high-order (HO) Blatter/Pattyn model (Blatter, 1995; Pattyn, 2003), a first-order approximation of the full-Stokes equations that assumes hydrostatic balance. The model incorporates stress-balance and mass-transport modules (Larour et al., 2012) and 3D thermal capabilities, that include advection (heat transferred by moving mass) and thermal diffusivity (Seroussi et al., 2013).

To adapt ISSM for Mars, we incorporate Martian gravity (3.71 m/s²), reasonable geothermal flux (0.025 W/m²; Parro et al., 2017), topography (Putri et al., 2019; Smith et al., 2001), surface temperature (150 K locked at the sublimation temperature of the overlying SPRC), and the basal topography mapped with SHARAD to give thickness (Putzig et al., 2018). MOLA 512 PPD data supplied the surface topography (Smith et al., 2001). In our latest simulations, for regions inside of the MOLA data gap, we use surface topography calculated with data from HRSC (Putri et al., 2019). Earlier simulations were completed before we became aware of this data set, so we performed various extrapolations to fill the data gaps (Table S1 in Supporting Information S1). The high latitude basal geometry cannot be constrained with any data set, so we extrapolate the base of the CO₂ deposits to have zero thickness at the extent of the boundaries made with geologic mapping (Phillips et al., 2011) and make assumptions about the thickness of the CO₂ units based on that geologic mapping (Table S1 in Supporting Information S1).

We also updated the material properties of the model fluid to match CO₂ ice (Table 1): rheology (Cross et al., 2020), heat capacity, thermal conductivity, and density (Piqueux et al., 2003). In these first-order simulations we did not include the top layer of water ice, the bounding layers (Alwarda & Smith, 2021) or the interpreted crevasses but anticipate their effects for future work (Figure S2 in Supporting Information S1). Once included, these bounding layers will likely reduce the modeled flow velocity.

The ISSM model domain for the ice cap is made up of a high resolution, anisotropic mesh of 13,500 finite horizontal elements and seven layers for a total of 94,500 individual elements. Mesh resolution cells range from 3.5 km (near the glacier margins) to 4 km (in the middle of the glaciers and in the no-data zones). The 3D thermal model is initialized with a conduction-only solution, forced by a surface temperature of 150 K and a geothermal heat flux of 25 or 10 mW/m². Based on laboratory experiments of CO₂ ice, thermal conductivity is set to a value of 0.4 W/m/K, and heat capacity is set to 700 J/kg/K (Ross & Kargel, 1998). We initialize ice rigidity (B) using the resulting 3D temperature field, according to a Nye formulation (Nye et al., 2000) updated with new experimentally determined constants for CO₂ (Cross et al., 2020), relating temperature to the flow law parameter A, where B = A^−1/n and has units of s^1/8 Pa (since n is set to a value of 8 for CO₂) and the activation energy of 67 kJ/mol. We disallow sliding (setting the friction coefficient to 500 [m/s]^−1/2) and set the initial thickness everywhere to 1 m (smallest permitted in the model).

To test sensitivity of the model to a variety of factors, we ran simulations with several variables. One variable we adjusted is a geothermal heat flux of 25 or 10 mW/m² fixed throughout the simulation. This allows us to test for the invert-ability of flow rate to find a constraint on geothermal flux. To test flow rates, we completed simulations in which flow was disabled or the stress exponent was varied (n = 8 or n = 6). n = 6 was chosen to split the difference between two solutions by Durham et al. (1999) with n = 5.6 (“recommended rheology”) and n = 7 (“plausible best estimate”). We also varied the domain size to test how the initial distribution of accumulation affected the final distribution of ice. See Table S1 in Supporting Information S1 for a full list of differences between the model runs.

Our deposition area for Simulation 08–11 is approximately 17:8 larger than that modeled by Buhler et al. (2019), and so we divide the surface mass balance (SMB) by that ratio to match the total CO₂ budget for the planet. We follow this ratio rule for the other simulations in which we varied the domain size. This implicitly assumes that CO₂ is not depositing in large quantities elsewhere, justified because the presently observed albedo and emissivity of the SPRC are significantly higher than other (seasonal) CO₂ deposits on Mars (Paige & Ingersoll, 1985), belying their locations and stability.

In two simulations, we modeled the current topography and ice thickness to estimate the speed of glaciers today. The geothermal flux was set to 10 and 25 mW/m² for simulations 12 and 36, respectively. Unlike the historical simulations, these were only run for 1,000 years to establish model convergence.

Finally, for each setup we solve a full transient solution (forward run) that includes a thermal solution with convection plus advection. The only varying input during any simulation is the SMB from Buhler et al. (2019). The transient solves the mass transport, stress balance, and thermal models through time, resulting in 3D solutions fields for ice thickness, ice velocity, and ice temperature. Inputs to the model are in Table 1 and Table S1 of Supporting Information S1, and results are discussed in the next section.

5 Glacial Model Results

During our model runs, CO₂ ice was deposited uniformly over the entire model domain, and we measured the total thickness of the ice in each cell as well as the temperature and three-component velocity. We also tracked the global ice mass and the maximum values of thickness, temperature, and velocity globally.

In order to have a basis for comparison, in Simulation 19, we ran the model with flow turned off and found that the entire model domain accumulated and lost ice uniformly, as expected. Total ice volume and ice volume per horizontal grid cell tracked the SMB, either increasing or decreasing, when SMB was positive or negative. The final result was a uniformly 44 m thick deposit with no lateral variation (see Supporting Information S1 and Movies S1–S6). This result did not match the observed distribution of MCID in any way.

For the models with flow turned on, no matter the variable being tested, the ice flowed downhill following the gravitational stresses and topographic gradient. Due to the large stress exponent governing the stress response of CO₂ ice (n = 8 compared to n = 3 for H₂O ice; Cross et al., 2020), ice that is deposited on steeper slopes experiences much faster flow than on lower slopes, even if it is thinner or warmer. We compared Simulation 20–24 with n = 8 to the “no flow” Simulation 19, which was otherwise identical, and found that the icy distribution was much closer to that observed by Putzig et al. (2018) and Alwarda and Smith (2021).

5.1 Stress Exponent

To test the sensitivity of the model to stress exponent, Simulation 25–26 was run with a stress exponent of n = 6 (a possible solution from Durham et al. [1999]) but otherwise identical to Simulations 20–24 and 19. The horizontal velocity near steep slopes dropped by nearly 30% from a max of 0.67 m/s to a max of 0.48 m/s. Likewise, for maximum thickness, the simulation with n = 6 had a much lower value of ∼360 m compared to the thicker ∼560 m for n = 8. The total volume at the end of the simulation was 9% lower for the n = 6 simulation than for the otherwise equivalent n = 8 case. Because the results of the higher stress exponent simulations were a closer match to the observations of Mars, we interpret this to support the laboratory measurements of Cross et al. (2020).

From all of our simulations with flow turned on, we find that the high stress exponent acts to focus the CO₂ ice into pre-existing basins where it reached more than 500 m thickness except in the n = 6 case (Figures 5 and 6, Movies S1–S6), leaving many locations with thinner ice. We interpret this process to be analogous to a lake serving as a catchment for rainfall. Rain accumulates only a few cm over a large area, but the level of the lake that catches this entire volume has a smaller area and will experience a much greater rise (Figure 7d). Compared to flowing water ice on Earth, the ∼63% reduction in gravity on Mars partially offsets the higher density of CO₂ ice (∼1,600 vs. 980 kg/m³ for water ice), so the stresses per unit thickness are 60% lower than what H₂O glaciers experience on earth. However, the weaker rheology of the CO₂ ice on slopes due to the higher stress exponent more than compensates for the lower gravitational stresses, and flow velocities reach several tens of cm per year for the steepest slopes during times of greatest SMB (Figures 5 and 6, Supporting Information S1). Locations with low slopes, including where ice is present today, had much more subdued velocities.

For all simulations with flow turned on, the greatest flow rates are always near steep slopes. Additionally, they always occur during periods of positive SMB, lagging slightly behind periods of maximum positive SMB because accumulation is ongoing (Figure 6). The maximum velocity reached for any point over the entire domain and duration exceeds 0.65 m/year. Velocity may approach 0 m yr⁻¹ for multiple reasons: locally for periods in which the ice is very thin and SMB is negative, and on low slopes.

6 Discussion

Previous work suggested that cycles of CO₂ deposition and sublimation built up the alternating layers of CO₂ and H₂O ice that compose the MCID (Bierson et al., 2016; Buhler et al., 2019; Manning et al., 2019). Bierson et al. (2016) and Manning et al. (2019) proposed that H₂O bounding layers form during periods of CO₂ sublimation, eventually cutting off the sublimation to preserve the MCID, but they did not propose how the bounding layers were deposited. Bierson et al. (2016) also observed that the CO₂ ice was preferentially thicker in depressions, but did not provide a physical mechanism for this observation. Buhler et al. (2019) showed that sublimation lags form H₂O bounding layers during periods of CO₂ ablation. That study was able to match the thickness ratio of MCID layers at one specific location, but made no attempt to explain the 2D or 3D distribution of the CO₂ ice. Neither model explained why the CO₂ ice concentrated in the valleys, how the final volumetric distribution occurred, or the unit thickness and stratification outside of one location. Alwarda and Smith (2021) found stratigraphic evidence supporting the Buhler model in the sequence of bounding layers but didn't offer any further explanations for the mass distribution.

Our modeling addresses all of these open questions. In particular, our results explain the spatial and volumetric distribution of the MCID that was heretofore unexplained. While depositing, the ice flows downslope into the depressions (Figure 7), creating the thickest deposits that can survive periods of negative SMB. The thinned deposits that lost mass to flowing ice do not survive sublimation periods, leaving the basement material exposed. Based on our Simulation 19, with flow turned off, the total SMB over 600 kyr is positive, but only enough to grow a few tens of meters of ice in total, uniformly over the domain, far short of the 400–1,000 m thickness estimated for the thickest observed deposits, and incompatible with the zero observed thickness in many locations within the accumulation domain. Thus, by focusing the ice into depressions, glacial flow explains the regions with thinner and thicker distributions and the persistance through periods with negative SMB, something a model without flow or highly preferential deposition cannot do. In that vein, if the flow speed is effectively reduced (via lowering the stress exponent or lowering geothermal flux) then less ice ponds to survive the negative SMB periods, and the total volume of ice at the end of the simulation is reduced. However, we found that changing the geothermal flux has a smaller effect than anticipated, likely because the thicker portions that have warmer basal temperatures are in basins (with shallower slopes), and the thinner, steeper sloped regions don't experience strong temperature increases that would accelerate flow.

Our modeling results demonstrate that glacial flow focused wide-scale accumulation into the curvilinear basins during periods of positive SMB and that the massive ice persists for long durations only in locations where the ice accumulates to thicknesses greater than the amount removed during the subsequent sublimation period with negative SMB (Figure 7c).

Simulations with a wider accumulation domain result in a worse fit between modeled and observed mass distribution, namely, the resulting MCID is thinner in important regions and more widespread, filling basins that have no CO₂ ice present today. This represents a major divergence that goes away with smaller model domains and supports the hypothesis of Colaprete et al. (2005) who found that atmospheric dynamics created the SPRC general distribution that we see today, a likely indicator of where the MCID deposition occurred over the past 600 kyr. Second, this supports the hypothesis from Buhler et al. (2019) that the SPRC and MCID are linked in formation and evolution (and spatially linked) because a wider accumulation of CO₂ ice over the duration of the model would have resulted in finding more MCID deposits that what has been observed, and thus a more widespread SPRC.

7 Conclusion

In this work we have identified numerous geomorphic features that are indicative of glacial flow of CO₂ units at Mars' south pole. The gross morphology and topographic contours of the CO₂ deposits, plus smaller features of crevasses and compression ridges, resemble the morphology of terrestrial glaciers.

To support the glacial interpretation, we updated the Ice-sheet and Sea-level System Model to work with martian conditions and flowing CO₂ ice using newly updated rheological laws (Cross et al., 2020). Using these setup conditions, along with known topography and unit thickness, we input the SMB from the most recent 600 kyr as modeled by Buhler et al. (2019). We find atmospherically deposited CO₂ ice will flow into existing topographic basins where it will pond. The final distribution generally corresponds in location and thickness to observed deposits (Putzig et al., 2018).

We ran models of varying distributions of ice deposition and ice flow laws and found that (a) simulations with reduced flow rates cannot reproduce the observed distribution of CO₂ ice and that (b) if the ice is deposited over a larger area, then the final result diverges in thickness and distribution from what we observe. The results support those of Buhler et al. (2019) who modeled the SMB over time and of Colaprete et al. (2005) who found evidence that supports a domain of deposition consistent with the current SPRC.

In summary, glacial flow of atmospherically deposited CO₂ ice explains the observed volumetric distribution of the MCID at the south pole of Mars, and glacial flow has been responsible for focusing widespread CO₂ accumulation into the curvilinear basins, where it has persisted for ∼510 kyr, even through periods of net loss. Model results suggest that the deposits are flowing presently but at a ∼1/3 reduced rate from their high flow periods, when SMB was at its maximum, at ∼4.8 × 10⁵ years before present (Figure 6). Glacial flow thus provides context for the resulting location, thickness, surface properties, and flow rate of the observed CO₂ deposits along with providing supporting evidence as to the origin and position of the SPRC.

Other exotic glacial landforms have been found in the Solar System, with the detection of an N₂ ice sheet near the equator of Pluto (Stern et al., 2015). This feature appears to have a very young surface age and exhibits evidence for thermally or sublimation driven convection cells (McKinnon et al., 2016; Morison et al., 2021). Other bodies are suspected to have thick, plastic, exotic ice deposits on their surfaces such as the moon Umbriel (Sori et al., 2017). With the increasing evidence for bodies of exotic ice present in the Solar System, there is a greater need to perform laboratory measurements to determine physical constants of exotic ices and to develop glacial models capable of quantitatively testing them. This in turn will test observations and interpretations, and enable us to invert for their climatic histories.