The Nature and Origin of Charon’s Smooth Plains

that the disclaimers , Abstract Charon displays extensive plains that cover the equatorial area and south to the terminator on the sub-Pluto hemisphere observed by New Horizons. We hypothesize that these plains are a result of Charon’s global extension and early subsurface ocean yielding a large cryoﬂow that completely resurfaced this area leaving the plains and other features that we observe today. The cryoﬂow con-sisted of ammonia-rich material, and could have resurfaced this area either by cryovolcanic eﬀusion similar to lunar maria emplacement or a mechanism similar to magmatic stoping where lithospheric blocks foundered. Geological observations, modeling of possible ﬂow rheology, and an analysis of rille orientations support these hypotheses.


Introduction
New Horizons has provided the first geologically useful images of Pluto's moon Charon (Stern et al., 2015), which have revealed that its geology is strikingly diverse . The Pluto-facing hemisphere, which was seen in the most detail by New Horizons during its flyby, is divided into two primary A C C E P T E D M A N U S C R I P T regions: Oz Terra and Vulcan Planitia. 1 Oz Terra consists of tectonically disrupted terrain that extends to the north pole and beyond onto the anti-Pluto hemisphere . A prominent series of scarps near the equator separates Oz Terra from a topographically low plains unit, Vulcan Planitia, that occupies at least 400,000 km 2 on the visible encounter hemisphere, and extends an unknown distance beyond the terminator (Figure 1 and units Sm 1 , Mt, and Rt in Figure 2). Vulcan Planitia contains craters, rilles, mountains, and other topographic features (Figures 2, 3, 4, and Robbins et al., 2018). This region exhibits a variety of textures: a smooth, lobate texture, a ridged/lineated texture, 1 Some names on Charon are now formalized and others are still informal, and follow the informal names in Moore et al. (2016) and Beyer et al. (2017). Vulcan is now referred to as a planitia instead of a planum after elevation data in Schenk et al. (2018).  and mottled textures.

A C C E P T E D M A N U S C R I P T
We present observations of these plains from the New Horizons Long-Range Reconnaissance Imager (LORRI, Cheng et al., 2008) and the Multi-spectral Visible Imaging Camera (MVIC) on the Ralph instrument (Reuter et al., 2008). This work analyzes these observations and examines possible analogs in the solar system, proposes possible mechanisms for a resurfacing cryoflow, estimates A C C E P T E D M A N U S C R I P T flow rheology, measures the orientation of tectonic features, and concludes that Vulcan Planitia is a unit emplaced by flow of icy material.

Observations
We define longitude and latitude on Charon according to the right hand rule and follow the recommendations of Zangari (2015). Charon's positive latitude pole points in the direction of the angular momentum vector and longitudes increase to the east. Charon's prime meridian is the sub-Pluto longitude. Pluto and Charon's poles are defined by Archinal et al. (2011a,b). Latitudes on Charon, following the current scheme for dwarf planets, are simply positive or negative.
Charon exhibits a visually prominent belt of fractures and scarps that define the boundary between Oz Terra and Vulcan Planitia. These fractures are highlighted by the solar illumination angle at the middle latitudes where they occur. The morphologic features in Vulcan Planitia are also evident due to the lower solar illumination angle. There is no information on structures southward of roughly −35 • due to the inclination of the Pluto system to its orbit around the Sun. All of that territory was in darkness during the New Horizons flyby.
New Horizons obtained extensive stereo coverage of the encounter hemisphere, and those data were processed to derive elevation data from the images (Figure 4 and Schenk et al., 2018) upon which our topographic interpretations are based. This terrain model has a ground scale of 300 m/pixel, and stereo height accuracy from ∼1.4 km to as good as 100 m in some places.
The mean elevation of Vulcan Planitia is one to two kilometers below that of Oz Terra, and the amplitude of topographic variations is much smaller across the plains than in Oz Terra . There are a host of possible explanations for these topographic amplitude differences, which could include regolith infilling the lows, or viscous relaxation affecting the landscape. However, little evidence has been identified for these two hypotheses. While there are landslides that have been identified within Serenity Chasma (Beddingfield et al., 2018, and unit La 1 in Figure 2 and Robbins et al., 2018), there are no observable geomorphic features that indicate infilling or movement of substantial additional amounts of loose regolith on Vulcan Planitia. If viscous relaxation were responsible, then we would expect other viscous relaxation signatures (such as up-bowed crater floors) that would be present in Vulcan and not in Oz, especially in similarly-sized craters in the two regions, and we do not. Therefore, we hypothesize that the lower topographic amplitude found on Vulcan Planitia is due to the plains-forming material itself infilling previously existing topographic lows. Supporting this interpretation, the lack of high-standing topography probably implies extensive subsidence of the original ice lithosphere prior to, and perhaps associated with, emplacement of the plains material.
Most of the highest-resolution images presented here were obtained with the New Horizons LORRI camera. The LORRI point spread function has a full width at half maximum (FWHM) of ∼2 pixels, which meets its specifications (Cheng et al., 2008), but can cause noticeable blurring. This blurring can be A C C E P T E D M A N U S C R I P T reduced through Lucy (1974) -Richardson (1972 deconvolution, and the images that this has been applied to are noted in the figure captions.

Surface Texture
The surface of Vulcan Planitia is smooth on 100-km scales , but its surface is marked with many smaller-scale features ( Figure 2). Numerous low hills, several kilometers across and several hundred meters high, occur on the planitia, often in clusters ( Figure 5). There are also broader domes and swells that produce a gently undulating surface over most of the planitia. In several places, particularly near the southernmost visible part of the planitia, some of these swells form convex sinuous features ('broad warp' features in Figure 2) associated with depressions ( Figure 6 and units Dm 1 in Figure 2): these convex margins appear similar to the convex margins of the montes (units Mo b and Dm 2 in figure 2 and detailed below in Section 2.3) and may similarly represent converging flow fronts as suggested by Schenk et al. (2018).
Some parts of the planitia are densely covered in mottled terrain with a locally-varying characteristic wavelength of about a kilometer ( Figure 5 and unit Mt in Figure 2). In some places there are sharp boundaries between these A C C E P T E D M A N U S C R I P T mottled plains and nearby unmottled plains, although the nature of the contact is uncertain .

Rilles
Vulcan Planitia is criss-crossed with a pattern of rilles, identified as small grooves (Stern et al., 2015;Moore et al., 2016 and the 'groove' feature in Figure  2 and in Robbins et al., 2018), which range from less than a kilometer to 3 km in width (Figure 7). These rilles range from a few kilometers to hundreds of kilometers in length, and in some locations show an arcuate or sinuous pattern. The resolutions of available images are insufficient to determine rille cross-sectional profiles in detail, although some appear to be simple graben, due to their inward facing walls. of meters deep. In some locations, these features display an en echelon pattern and exhibit Y-shaped intersections, consistent with a graben interpretation. These rilles appear to cut, and thus post-date, many of the broad domes and swells. However, the available images do not display unambiguous examples of grooves cutting any of the numerous impact craters that are superposed on the plains unit, or their ejecta. Therefore, the grooves likely formed soon after plains emplacement, before subsequent cratering events. Rilles and mottled terrain of Vulcan Planitia visible in this mosaic of LOR 0299180418 and LOR 0299180421, from the C MVIC LORRI CA observation, deconvolved to enhance features. The image is centered at approximately 358 • E, +3 • . Stern et al. (2015) and Moore et al. (2016) described Vulcan Planitia's mountains. There are two large solitary massifs (Butler and Kubrick Mons) and one cluster of smaller peaks (Clarke Montes). Vulcan Planitia does not cleanly onlap the toe of these mountains, but instead displays a convex, depressed margin that surrounds them (e.g. Figures 8 and 9), engendering the term 'mountain in a moat' for these features.

Mountains
Butler Mons is the largest massif ( Figure 10), but by the time of closest approach it had set below the terminator. Butler can be seen in New Horizons approach images with substantial scarps to its east and south. These scarps may indicate a hitherto unknown southern boundary to Vulcan Planitia that 8

ACCEPTED MANUSCRIPT
A C C E P T E D M A N U S C R I P T might mirror the Oz Terra scarps to its north. Kubrick Mons is also a single massif that measures 30 by 20 km in plan form, and rises 3 km above the Vulcan Planitia surface with a 1 km deep encircling moat (Figures 8 and 9). Kubrick displays a ridge crest and other sharp-looking ridges that run down its flanks. This morphology on Kubrick (and to a lesser extent on Clarke) is similar to morphologies observed in Oz Terra, and not similar to the surrounding rounded lobate shapes of Vulcan Planitia. This morphology implies that Kubrick (and other Mo b units on Figure 2) is older than the Vulcan Planitia surface.
There are also at least two other, smaller features that share a strong resemblance to Kubrick. They have steep sides, and a depressed margin ( Figure 11 and other unit Mo b outcrops in Figure 2).
The Clarke Montes is a cluster of peaks that sit together in a single, variable margin depression ( Figure 12). They stand ∼1 km above the surrounding plains. The northern margin exhibits a 2 km drop from Vulcan Planitia, but the southwestern margin exhibits only a 1 km drop. The east margin does not have a distinct lobate margin, but appears to smoothly grade out from the floor of the depression up to the plains ( Figure 13).
Two hypotheses have been proposed for the moats which surround these montes : (1) downward flexure of pre-existing plains material from loading by younger mountain material (formed for example by extrusive cryoflow), or (2) incomplete embayment of pre-existing mountains by the plains material.
The downward flexure hypothesis was explored by Desch and Neveu (2017) and they came to a conclusion similar to one discussed internally by the New Horizons team that the thickness of the elastic lithosphere, T e , would need A C C E P T E D M A N U S C R I P T Figure 10: These LORRI observations of Butler Mons show it before it sets beyond the terminator, but its peak remains illuminated in (e).
Shadows indicate a height of 4.5 km .
The arrow is offset and points to the peak of Butler Mons in all figures.
All of these are LORRI observations: (a) PC MULTI MAP B 18 02, (b) PCNH MULTI LONG 1D1 02, (c) PC MULTI LONG 1D2A 02, (d) C LORRI FULLFRAME 1, (e) C LORRI . These are not map projected, but north is generally to the bottom of the frame and east to the left.  to be very thin (less than 2.3 km). Desch and Neveu (2017) hypothesized that the warm icy material that erupted onto the surface was from a local hotspot that was warming the subsurface in this region, and thinning the elastic lithosphere in the area near Kubrick. The Desch and Neveu (2017) hypothesis is consistent with the rift flank T e estimates obtained by Beyer et al. (2017), but we concluded such a low T e was unlikely to be correct given the fact that large craters do not display evidence of viscoelastic relaxation, e.g., up-bowed floors. Additionally, the morphology of Kubrick doesn't resemble what might be expected of a constructional cryovolcano, for example, central depressions consistent with a caldera, etc. The embayment hypothesis is supported by the presence of similar convex scarps elsewhere on Vulcan Planitia surrounding depressions without mountains, as noted by Moore et al. (2016). Furthermore, the arrangement of Clarke Montes and the presence of a break in slope at the base of its enclosing scarp, which is in some locations several kilometers from the massifs, would be consistent if the scarp was a flow front, but not a flexural bending of the lithosphere. It is not clear, however, why the flow would have stopped close to the mountains but not at their base. Perhaps the base of the Clarke Montes changes from the steeper slope of the montes and shallows, but still slopes away from the mountains, causing an oncoming flow to stall before reaching the mountains themselves. Also, the local hotspot hypothesis does not explain the depressed margin around the Vulcan Planitia boundary, but the embayment hypothesis does.

Margins
Along the south-facing boundary scarps of Oz Terra  that define Vulcan Planitia's northern boundary, the surface of the planitia is depressed by a kilometer or more in elevation within about 20 km distance of the boundary (Figures 4, 14, and 15). Where Vulcan Planitia meets the scarps of Oz Terra, the Vulcan Planitia units appear to onlap the scarps, consistent with a surface that was emplaced after the scarps formed.
In the portion of Vulcan Planitia between Spock crater and where Serenity Chasma's eastern end transitions from a graben to a south-facing scarp, there is a well-resolved pattern of closely spaced and parallel rilles and ridges, which parallel the boundary. This surface texture only appears at distances greater than 15 km from the bounding scarp where the planitia materials are presumably thicker, implying a pile-up of flow material where it was thick enough to do so. Borrelli and Collins (2018) examined the shapes of the margin depressions along the Oz/Vulcan boundary, and tested a variety of viscous flow and plate bending models. Their calculations did not strongly favor any mechanism, nor were any disallowed. They conclude that the morphological evidence better supports the advancement of a spreading viscous cryoflow instead of a plate bending model.

Analysis and Results
Charon is comparable in size and density to several icy satellites of the giant planets, some of which also show resurfacing in the form of plains (although to different extents). Figure 16 compares the masses and densities of these satellites, indicating whether or not resurfacing is observed on their surfaces. The approximate location of the toe of the south-facing slope into Vulcan Planitia, estimated from the images, is shown by the end of the shading. Just to the south (right) of that point, there is an area that is depressed by about a kilometer for a length of about 20 km, before it rises up to the level of the larger Vulcan Planitia. The average elevation of Vulcan Planitia can be seen to increase, and the depth of the margin decreases, as one goes eastwards. In this area of the terrain model, the stereo height accuracy is a few hundred meters, and wouldn't plot as a visually separate line on these graphs.
In general, higher satellite densities and masses appear to be associated with resurfacing, although not uniquely.

Rheology of the Flow
We interpret the lobate, rounded features on Vulcan Planitia to be the result of flow emplacement of the surface. The lobate margins, which encircle the mountains, appear to be the result of viscous flow that has encountered an  Figure 16: This plot of density versus mass for the large icy satellites in the outer solar system shows that Charon is typical among this group. Solid symbols indicate those objects that have some evidence of resurfacing, question marks indicate that there is not enough data to determine whether there is any resurfacing. Source data and references for this figure can be found in Table 4.
obstacle. The depressed margins along the Oz Terra scarps also imply a viscous flow, but the difference in width between the moats around the mountains and the depressed margin along the scarps may be due to the different depths of the viscous fluid in the center of Vulcan Planitia versus the potentially shallower depths at the margins.
To estimate a viscosity for such a flow in Vulcan Planitia (although we envision many happening simultaneously), we will make a number of simplifying assumptions: that the flow was emplaced from a single eruption event, that the viscosity was high enough that inertia can be neglected, and that the pre-emplacement surface was horizontal (although we suspect that it was not). Under these assumptions, the emplaced fluid spreads under its own weight and the equations of Huppert (1982) can be used to analyze its motion.
For a two-dimensional flow, the dynamic viscosity, η, can be derived using Here η N is a dimensionless constant which is equal to 1.411 for a flow that is emplaced instantaneously (Pattle, 1959), ρ is the density, h is the flow thickness, L its length, g the acceleration due to gravity, and t the time that the flow was in motion. The approximation results from taking the cross-sectional area of the flow to be hL. Thinner or faster-moving flows imply lower viscosities.
The quantity η N changes slowly as the flow emplacement condition is varied (Huppert, 1982), so even taking this quantity to the fifth power contributes only a minor amount of uncertainty. The greatest uncertainty associated with equation 1 is determining t, the duration of flow. The duration of flow is a function of the rate that the material cools, and therefore the time that it takes the material to freeze. Schenk (1991) suggests two possible time scales for the rate of cooling. The first is the conductive timescale t cond : where κ is the thermal diffusivity. This timescale is appropriate if the material has a high enough viscosity that vertical mixing is minimal. The second is the characteristic radiative cooling time, based upon a lava cooling model by Crisp and Baloga (1990) and is appropriate if the interior of the flow is wellmixed and radiation from the warm interior (e.g. via cracks) dominates over conduction. The (Crisp and Baloga, 1990) model assumes fractional areas of exposed material between 0.001 and 0.1 based on values observed for several types of terrestrial lava flows. In this case the timescale is where C p is the specific heat capacity, is the emissivity (taken to be 1), σ is the Stefan-Boltzmann constant, and T 0 is the characteristic interior temperature. The ratio of these two timescales is then T 3 0 σh/k, where k is the thermal conductivity. On the Earth, the conductive timescale greatly overestimates observed flow durations, while the radiative timescale is more accurate (Schenk, 1991). The discrepancy between the two timescales, however, will be significantly smaller on Charon than on Earth, because the likely value of T 0 is so much smaller. Table 1 gives thermophysical properties for four potential cryo-fluids.
As an example, assume a flow with h = 2 km, L = 30 km, and g = 0.288 ms −2 . The lobate features interpreted as flow fronts along the Vulcan

A C C E P T E D M A N U S C R I P T
Planitia margin and the moats around the montes have a 1 km elevation, we assume 2 km here to account for a slightly deeper pre-Vulcan surface. We assume a 30 km length because we don't see morphologic evidence for long distinct flows, but this parameter could assume a number of values. Here we take T 0 to be set by the melting point of the substance, T melt , and the results are given in Table 2.
The viscosities derived assuming radiative cooling (∼10 13 -10 15 Pa s) are much higher than likely liquid or slush viscosities, but reasonable for glacial ice on the Earth (which is at the low end of that range) and viscosities derived by Jankowski and Squyres (1988) for flows on Ariel (at the upper end of that range). The viscosities derived assuming conductive cooling (∼10 16 -10 17 Pa s) are similar to those derived for solid creeping water-ice, but unreasonably high for N 2 and CO so close to their melting points (which rules them out as likely cryo-materials in this scenario).
Due to the much higher latent heat of fusion for likely cryoflow material candidates vs typical silicate materials, κ in equation 2 might need to be replaced by a value that is effectively 10 to 20 times less than κ, which would increase the values in Table 2 by an order of magnitude (O. Umurhan, pers. comm.).
In any case, such high viscosities suggest that very little vertical mixing will occur, and thus that conductive cooling is an internally self-consistent assumption. Whether other assumptions (e.g. emplacement as a single flow) are correct is currently unclear.
If we assume that Kubrick and Clarke are surrounded by Bingham rheology cryoflows, and that the Vulcan Planitia margin depression is also the result of a Bingham rheology cryoflow that was flowing in a direction from the center of Vulcan Planitia outwards towards the bounding scarps that define the border with Oz Terra, we can use the profiles in Figures 9, 13, and 15 with the following equation, modified from Melosh (2011, equation 5.11): to estimate the Bingham yield strength (Y B ). The thickness of the flow is H, and x is the distance from where H is measured to the edge of the flow. When we do, we find yield strengths from 8.3 × 10 3 to 4.6 × 10 4 Pa. These are similar to the yield strengths reported by Melosh and Janes (1989) for flows on Ariel, and are in the range of yield strengths for terrestrial silicate lavas (Moore et al., 1978).
There are a variety of other, more complete models of lava flow rheology that could be employed (e.g. Sakimoto and Gregg, 2001), but our data do not resolve flow margins or levees, nor do they reveal other details that would be important input to models of that nature. Additionally, the results presented here make it difficult to distinguish between solid-state extruded flows and a viscous mixture of liquid and crystals (Kargel et al., 1991;Schenk, 1991), especially in view of the rather poorly-characterized properties of low temperature, multi-phase mixtures (Zhong et al., 2009).

Rille Tectonics
The rilles within Vulcan Planitia formed either contemporaneously with or following plains emplacement. The morphology of the rilles in Vulcan Planitia is consistent with a tectonic origin, due to their regular spacing, mostly linear character, and negative topography. In addition, the en echelon arrangement of these rilles are consistent with extensional graben . If the rilles are extensional graben, then their orientations are determined by the stress patterns present during their formation. Consistent rille orientations across Vulcan Planitia would reflect a consistent and large-scale stress field. Charon likely underwent extensional stresses following its formation (e.g. Beyer et al., 2017).
Rilles present within Vulcan Planitia appear to exhibit consistent orientations in adjacent regions. However, additional rilles, with more complex orientations, are present in regions farther from the Oz Terra boundary. Here, we compare rille orientations across Vulcan Planitia to quantify their spatial variations in azimuth (see the appendix for more information). To conduct this analysis, we measure and compare the distribution of rille orientations within various regions (Figure 17).
Nearly all analyzed areas exhibit prominent rille azimuths that parallel the orientation of Serenity Chasma. Eleven out of the 14 analyzed areas exhibit prominent E-W rille orientations, and this orientation is the most prominent across Vulcan Planitia. Many areas exhibit rille azimuth distributions similar to their proximal neighbors ( Figure 17). For example, Areas 4, 6, and 8 exhibit rilles more similar in distribution to each other than with those near Serenity Chasma, with a prominent N-S orientation. The variation in Vulcan Planitia rille orientations with distance from Serenity Chasma suggests that another stress mechanism, in addition to global expansion, has acted on the lithosphere in these regions during rille formation. It is worth noting that there seems to be an absence of rilles in the areas surrounding Kubrick and Clarke Montes (also noted by Robbins et al., 2018). This lack of rilles may be the consequence of local stresses associated with these peaks, or from stresses related to the resurfacing cryoflow during emplacement or cooling.
We estimate the total strain contributed to Vulcan Planitia from rille formation. Because the rille widths are too small to accurately measure using available images, the strain for each rille was estimated based on their lengths. The maximum displacement-to-length ratio for extensional faults is 0.02, and A C C E P T E D M A N U S C R I P T Figure 17: Rose diagrams showing the spatial variation in Vulcan Planitia rille orientations. Each rose diagram is normalized, so that the first rille azimuth mode is represented as the longest sector, and the lengths of the remaining sectors are scaled proportionally. The N values on each diagram represent the number of rille segments in each area. The most prominent rille orientation is E-W, which parallels normal fault scarps within Serenity Chasma to the north. Other prominent rille azimuths, consistent between adjacent areas farther from the Vulcan Planitia boundary, are NE-SW and N-S. This data is in a simple cylindrical projection, from +30 • to -30 • latitude and 300 • E across the prime meridian to 50 • E. For more information about our methodology, see the section in the appendix starting on page 33, and supporting data for the rose diagrams can be found in Table 5.

A C C E P T E D M A N U S C R I P T
this ratio has been measured on icy satellites, including Europa (Nimmo and Schenk, 2006) and Charon . The increase in surface area for each rille is the sum of the two inward facing exposed fault scarp areas that make up the graben walls, scaled by the cosine of the dip angle. The average maximum rille displacement is ∼ 1 km, assuming the rille walls exhibit 60 • dips. The total exposed new rille planform area is estimated to be 6,680 km 2 . For an estimated Vulcan Planitia area of 400,000 km 2 to 500,000 km 2 , this areal expansion corresponds to an estimated horizontal strain range of 1.4% to 1.7%.

Resurfacing
Our observations and analysis above indicate that the area of Vulcan Planitia was resurfaced as compared to the more 'primary' surface in Oz Terra. However, the manner of this resurfacing is difficult to ascertain. We present three hypotheses below: (1) outright melting and refreezing of the lithosphere in place (Section 3.3.1), (2) lunar maria-style flood volcanic cryoflows (Section 3.3.2), and (3) foundering of lithospheric blocks into a viscous mantle (Section 3.3.3). We favor the latter two hypotheses because the presence of embayed or submerged, intact blocks supports some of the observed features that melting does not, and the melting-in-place hypothesis has some thermodynamic and energetic challenges that would need to be resolved. The cryoflow (Section 3.1) and rille (Section 3.2) investigations discussed above are independent of these mechanisms, and mostly explore what these surface features can tell us about the character of the material that resurfaced Vulcan Planitia, regardless of what mechanism was responsible for causing it.
The tectonic structures in Oz Terra are fairly ancient as they cut few craters, yet are preserved via the lack of saturation impact cratering . The transition between Oz Terra and Vulcan Planitia appears abrupt in available New Horizons images, but we propose that these two terrains are natural terrain end-members. Between those two end-members are the tilted surfaces of the lithospheric blocks in Oz Terra that directly border Vulcan Planitia (Figure 4). Similar tilted blocks are observed on Ariel (Schenk, 1991;Schenk et al., 2018). These blocks on Charon represent a point on the continuum between Oz Terra lithospheric blocks that are not tilted, but translated away from each other, and the blocks that we think existed on the pre-Vulcan surface before the cryoflows.
The one to two kilometer difference in average elevation between Vulcan Planitia and Oz Terra  is difficult to explain. Malamud et al. (2017) without the benefit of knowing that Vulcan had a lower elevation proposed that perhaps Vulcan was a giant impact basin. That could be an explanation, but there doesn't seem to be any significant morphological evidence that it is the case. Alternately, one could imagine that there was a subsurface difference in properties that could have contributed to this variation. As the ice shell was freezing, it would be a thickening spherical shell, over a spherical ocean, over a cooling core (Bierson et al., 2018). If the poles were colder than the equatorial regions, perhaps due to insolation or early tidal dissipation energy still being deposited in Charon, this might cause the ocean to transition from

A C C E P T E D M A N U S C R I P T
spherical to a band beneath the surface as the last warm spots were centered on the equator. As the ocean froze, pressure could have built up in the remaining liquid (Manga and Wang, 2007). In this situation, when the ice shell cracked and fractured from the extension, the poles would have a strong 'support' but the equator would be weak, and lithospheric blocks would more easily sink or subside in that area. Furthermore, the pressurization may have helped push the remaining buoyant liquid upwards but once it erupted, the pressure was released and erupting material became harder to mobilize, leaving the Vulcan Planitia area just a little lower in elevation than Oz Terra. Unfortunately, there are currently no data that might suggest this situation. Insolation models of Charon based on current day inclination indicate that the opposite is true: that the poles receive more solar energy than the equator integrated over millions of years . However, those models don't reach into the first billion years of solar system history.

Melting In Place
The resurfacing of Vulcan Planitia could have been the result of a largescale thermal event centered on the equator in this region of Charon. This event could have caused the existing lithosphere to melt in place at the time of crustal expansion. Such an hypothesis would indicate that Vulcan Planitia is simply the melted and refrozen remains of the pre-existing lithosphere.
In this scenario, the montes present on Vulcan Planitia could still be explained in one of two ways. They could be blocks that were not melted completely and simply bobbed and rotated into their present position and then were frozen in place, implying a mostly water-ice composition such that the lithospheric blocks float on the melted material. Alternately in this scenario, their explanation as constructional features over local hotspots might be strengthened and explained as these hotspots being the last remnants of this large thermal event.
The tilted blocks at the Oz boundary could be the expression of the edge of this thermal melting event, such that these blocks did not melt or sink, but with disruption at depth in this boundary zone sufficient to severely tilt them.
A similar hypothesis was put forward for chaos formation on Europa (O'Brien et al., 2002). It required a thin ice shell and tidal heating supplied by Europa's orbit around Jupiter. Charon is unlikely to have had a shell so thin at the time of resurfacing, and probably lacked the strong input of tidal energy. Finding a mechanism to provide a regional thermal pulse is challenging, and so is the amount of energy needed to melt so much near-surface material that is near 50 K today (or lower earlier in solar system history) and can radiate unimpeded to space.

Cryovolcanism
Voluminous eruptions of cryo-materials are to be expected on midsize icy satellites, due to the low driving stresses due to buoyancy (low gravity) and the presence of NH 3 (ammonia) hydrates or other molecules which could provide the needed density relationships between solid and liquid (Stevenson, 1982). This hypothesis advances that cryo-material was erupted onto the low, pre-Vulcan surface from the interior and cryoflows flooded the area, emplacing the Vulcan Planitia we observe today (Figure 18).

A C C E P T E D M A N U S C R I P T
The icy lithosphere of Charon broke into large polygonal blocks in response to global scale tensional stresses . We envision that this breakup was global and also extended to the area that was to become Vulcan Planitia. It is possible that this breakup was more severe in this area than in Oz Terra, or that whatever caused the pre-Vulcan area to be topographically

A C C E P T E D M A N U S C R I P T
lower promoted significant lithospheric weakening in this area that permitted cryo-material to rise up and erupt onto the surface.
The traditional difficulty to overcome with cryovolcanism is how to develop a density inversion on an icy world that would allow cryo-material to rise and erupt. The density properties of pure water ice would cause icy lithosphere to float on an icy mantle, even if it were physically disrupted. Johnson et al. (2017) suggests that salt content and porosity differences may provide such a density inversion. For Charon, we propose that the presence of an NH 3 hydrate mixture in Charon's lithosphere and mantle are the solution for getting a density inversion (Kargel et al., 1991;Schenk, 1991;Kargel, 1992;Hogenboom et al., 1997;Cook et al., 2007). It should be noted that the density differences between an NH 3 hydrate liquid and its solid (or even pure water ice) are not very large and an NH 3 -H 2 O peritectic liquid has a lower density than solid ammonia hydrate, but slightly greater than pure water ice (Hogenboom et al., 1997), but all of these densities are very close to one another, and small variations may allow extrusion of NH 3 hydrate even through a water ice layer.
The broad surface of Charon is primarily water ice, but Grundy et al. (2016) andDalle Ore et al. (2018) note that at least a few craters appear to excavate NH 3 hydrates from the subsurface. Cook et al. (2017) also find evidence for NH 3 on Pluto's satellites Nyx and Hydra. Based on these observations and considering the likely incorporation of NH 3 in Kuiper Belt objects (McKinnon et al., 2008), we have an expectation of a bulk abundance in the Pluto-Charon system of about 1% compared to water. However, ammonia would concentrate in any primary melt or residual ocean subject to freezing.
We hypothesize a possible case involving NH 3 hydrates that would achieve a density inversion and support the observations. We assume that beneath the ancient surface of Charon, a subsurface ocean formed and then refroze (Bierson et al., 2018). In this ocean, water would freeze out of the solution first, enriching the remaining 'ocean' in NH 3 as freezing progressed (Kargel, 1992). This process would leave a mostly original bulk composition (mostly water ice, but with rock and other contaminants) lithosphere overlying a water-NH 3 ocean, and would also concentrate NH 3 in that ocean, even from relatively small original bulk amounts. If the lithosphere were pure water ice, it would still float on a water-NH 3 ocean. However, as mentioned above, the densities are very similar, and if the lithosphere were not pure water ice but was an NH 3 hydrate that contained even a small fraction of rocky material, that might enable a density inversion. When the lithosphere expands and is disrupted, the water-NH 3 ocean material would rise up above the mostly-water-ice lithospheric blocks. This cryo-material wells up, resurfaces the area and then proceeds to freeze. This newly-emplaced unit should preserve the bulk composition of the 'ocean' that it came from, because it would have cooled much faster than the slow fractional crystallization described above. With such a high ammonia concentration, one might have expected the spectrometers on New Horizons to have detected it, but photolysis and proton bombardment act to remove this constituent from the surface over a few tens of millions of years (Cooper et al., 2003;Moore et al., 2007;Cassidy et al., 2010). Therefore, spectrometers on New Horizons

A C C E P T E D M A N U S C R I P T
would not have detected NH 3 hydrates in substantial amounts on the surface of Charon in this case, even though they are probably present not far below the surface. However, the fact that this spectral feature is seen in spectra of the ancient surfaces of Nix and Hydra Cook et al. (2017) casts some doubt on the anticipated destruction rate, or perhaps the identification of the feature with NH 3 or NH 3 hydrate.
It is important to note that there are different 'styles' of large igneous provinces in basaltic volcanism, and under this analogous hypothesis for Charon we are advocating a style that is similar to the volcanism of the lunar maria where there is an absence of clear vents or feeder dike systems, and it took very-high-resolution gravity data to reveal the volcanic plumbing of the maria (Andrews -Hanna et al., 2014). The other 'style' of large igneous province are the large sheet-flow flood basalts on the Earth like the Deccan Traps or the Columbia River Basalts. These mechanisms have a few discrete volcanic centers as sources, and have lavas that flow out over very large areas. Application of this style to Charon would require only a few discrete sources for the cryoflow, and then a large volume of cryo-material needs to stay mobile long enough to cover a large area, which seems thermodynamically challenging.
The lunar maria style allows for many 'local' sources of cryo-material all across the region, and also explains the lack of large, obvious flows. Moreover, the surface swells that form convex sinuous depressions (unit Dm 1 and the 'broad warp' features in Figure 2), could mark eruptive centers, if late-stage cryo-materials were withdrawn.
In this hypothesis, the montes on Vulcan Planitia are simply pre-existing areas of high elevation that were embayed by the cryoflows like kipukas or nunataks. These montes are not significantly different in elevation than the high terrain in Oz Terra and could simply be remnants of higher-standing topography in this area. The tilted blocks at the Oz boundary could be tilted due to mechanisms that acted in the pre-Vulcan area related to the subsidence of the pre-Vulcan surface and not related to the flows, or they could have been tilted due to subsidence when the cryo-material erupted from depth.
Charon need not have undergone extensive internal melting of its ices, but an ocean earlier in geologic time is plausible (Bierson et al., 2018), and it is a natural location for NH 3 to be sequestered. Even if non-buoyant compared with the porous icy lithosphere above, internal pressurization (Manga and Wang, 2007) may have led to the eruption of cryo-materials. We estimate that Charon should have accreted sufficient nitrogen (in the form of ammonia, or else converted to ammonia via subsequent chemical processing) to explain such an extensive cryovolcanic unit as Vulcan Planitia. For example, a 400,000 km 2 plains unit 2 km thick and of approximately dihydrate composition would require a bulk NH 3 abundance of ∼ 2 × 10 4 , about 20 times less than Charon's likely accreted NH 3 -ice (McKinnon et al., , 2008. Sequestration in a global residual ocean of dihydrate composition (32 wt% NH 3 ) would only yield an approximately 15 km thick basal layer, however, unless there is another source of NH 3 (i.e., core organics) or other substances (salts, methanol) are involved (McKinnon et al., 2008).

Foundered Crustal Blocks
The above hypothesis assumes a moderate disruption which fractures the pre-Vulcan landscape enough for cryo-materials to erupt. This new hypothesis is a modification of the previous hypothesis in which the disruption is so great, that the individual blocks of the pre-Vulcan surface lose their support and founder into the interior followed by the mantle material rising up under the same density inversion and chemistry arguments as above and resurfacing Vulcan Planitia (Figure 18). Kattenhorn and Prockter (2014) outline a subduction mechanism where large plates on the surface of Europa are subducted under one another driven by a plate-tectonic mechanism, supported by further work by Johnson et al. (2017). What we are proposing here is different from plate-tectonic subduction and is a foundering block mechanism (first proposed by Schenk et al., 2018) similar to terrestrial magmatic stoping. The lithosphere is extensively fractured, and the resulting blocks sink due to a density inversion and potentially rotate into the underlying material that wells upwards and envelopes those blocks. The erupted cryo-material flows across the surface of these submerged blocks and then freezes in place creating a new lithospheric surface. The isolated massifs identified above would in this case be the tips of such tilted polygonal lithospheric blocks, which stand above the plains of Vulcan Planitia like kipukas or nunataks, as above.
The tilted blocks on the border between Oz Terra and Vulcan Planitia under this hypothesis are explained as transitional between the Oz Terra blocks that are just translated away from each other, and the blocks that foundered, sunk, and were covered by Vulcan Planitia plains material. Of course, in this density inversion situation, when the solid is denser than the liquid, the resulting lithospheric forces will be compressive, and not extensional at the surface, which may seem to contradict the global expansion hypothesis, but this density inversion situation would only happen late in the freezing process, still allowing for global expansion, but may also point to incomplete ocean freezing.

Strain in the Plains
Irrespective of the exact mechanism, the emplacement of warm flows on the surface would result in subsequent cooling and contraction of those flows. As an example, a thin sheet of water ice cooling from 240 K to 40 K has a mean linear thermal expansion coefficient of about 2 × 10 −5 K −1 (Röttger et al., 1994). The areal strain resulting from the expansion of this cooling would then be about 0.8%. Water ice expands, but NH 3 hydrate will contract, and Croft et al. (1988) indicate that for ammonia water liquid, the volumetric thermal expansion coefficient is about 5 × 10 −4 K −1 . The areal strain resulting from the contraction of cooling from 240 K to freezing at 175 K is about 2%. This can be compared with the areal strain of ∼ 1.5% across Vulcan Planitia estimated in section 3.2. On the other hand, if the cryo-material flow is not a thin sheet, but is thicker, then there would also be vertical contraction of the emplaced material, perhaps up to several hundred meters. Given the substantial uncertainties in the properties of the emplaced material and in the inferred strain, this apparent agreement is quite encouraging.

Cryovolcanism Hypothesis
Foundered Block Hypothesis Subsurface ocean begins freezing, causing expansion of the lithosphere, and therefore breaks up the lithosphere. At the positive (north) pole, this expansion was moderate, and produced what is now Oz Terra (the negative pole might be similar, but we lack data). At the equator, the expansion was more severe and subsidence occurred, creating the pre-Vulcan lowlands.
At the equator, the disruption of lithospheric blocks was severe, and the blocks foundered and sank. The ammonia-rich mantle rose up through the disrupted lithosphere and resurfaced Vulcan Planitia. The cryo-material flowed out from many vents and covered the pre-Vulcan surface, embaying the existing highs and covering everything else. This process mostly buried the source vents, but the large depressions (units Dm 1 in Figure  2) could be vents where at the very end of the process, cryo-material was flowing backwards down into them, making them look more like sinks than sources.
The 'rising tide' of flowing cryo-material enveloped the sinking blocks to varying degrees. Some blocks completely sank, leaving no trace on the plains surface; some did not sink, and the cryoflow surrounded them, leaving the montes we observe. Some blocks sank just enough that the cryoflow surrounded and covered them, but as they sank, that flowing material froze solid before an equilibrium surface could be established, leaving the large swells observed. The boundary between Oz Terra and Vulcan Planitia is the result of whatever pre-Vulcan process created the elevation difference between Oz Terra and the pre-Vulcan surface. The tilted blocks on the boundary may have been tilted by that process or during Vulcan emplacement from the flow of mantle material from the subsurface.
The boundary between Oz Terra and Vulcan Planitia is the difference between lithospheric blocks that were covered (in Vulcan) and those that weren't (in Oz). The blocks on that margin that did sink, perhaps did so to a lesser degree than those in the 'center' of Vulcan Planitia.
At the Oz Terra/Vulcan Planitia boundary, the cryoflow would have been flowing out from the center of the Vulcan Planitia mass across either the subsided pre-Vulcan surface or across the tops of sunken lithospheric blocks towards the boundary. The flow would have encountered a shallower 'basin' near the boundary and this may be the reason for the depressed margin we see today. Part of that outward 'flow' resulted in the pile-up feature north of Spock (between profiles 9 and 11 on Figure 14). After the cryoflow had frozen in place, it would have undergone a cooling period during which Charon itself may still have been in the last stages of extension, and this cooling and extensional force may have caused the larger rilles in the plains (Hillier and Squyres, 1991) that parallel the Oz Terra boundary. After that, the only modification to the plains surface has been impact cratering.

Age and Sequence of Events
Cratering on Vulcan Planitia and other terrains on Charon is discussed in . This abundance of craters is consistent with an old age of ∼4 Ga found for Charon's northern and southern terrains Singer et al., 2016). This age estimate is based on crater distributions and modeled impact rates . It is notable that within Vulcan Planitia there are no unambiguous examples of fractures cutting craters, or other indication of endogenic modification or infilling of craters, suggesting that the resurfacing and further endogenic modification of Vulcan Planitia occurred rapidly relative to cratering timescales. Vulcan Planitia appears to embay Oz Terra from the south, and has at least one possible outlier within Oz Terra  and unit Sm 2 on Figure 2 and in Robbins et al., 2018). Vulcan Planitia is somewhat younger than the terrains of Oz Terra Robbins et al., 2017Robbins et al., , 2018, and we lay out a hypothesized relative timing of features and events based on the cryovolcanism and foundered lithospheric block hypotheses in Table 3.
That sequence assumes a single 'ocean freezing' in Charon's history to drive this activity. This is consistent with the possibility of a subsurface ocean (Bierson et al., 2018) or intense serpentinization (Malamud et al., 2017) resulting in a large expansion event at approximately 1 Ga after formation. Desch and Neveu (2017) have hypothesized that there might be two 'ocean' episodes in Charon's history. Their proposed long-lived initial ocean doesn't seem to agree with Vulcan Planitia formation early in Charon's history, which is constrained by the cratering record, but it does raise the possibility of different episodes in the expansion story above. The Malamud et al. (2017) modeling also indicates the possibility for a secondary, smaller period of expansion. It is possible that there was an initial episode that broke up Charon's lithosphere globally into large blocks. Then perhaps there was a significant pause after that initial disruption, and then a second episode (perhaps more localized to the equator or just this part of Charon) occurred and caused the cryo-material eruption or block foundering. Of course, these episodes can't be too far apart, or last too long, because Charon's cratering record shows that Oz Terra and Vulcan Planitia have very similar ages.

Conclusions
Large icy worlds in the solar system exhibit resurfacing to different degrees. Other objects for which resurfacing is observed either have limited areas of resurfacing, later activity that overprints the resurfacing, or just poor data resolution, such that developing and then testing hypotheses are difficult. We argue that the global expansion of Charon's surface  and subsequent resurfacing were both the result of the progressive freezing of a subsurface ocean. This provided a surface that stands as a witness plate to this geologic resurfacing, without much overprinting by other processes, which was excellently observed by New Horizons (Stern et al., 2015).

A C C E P T E D M A N U S C R I P T
We hypothesize that the resurfacing of the plains is not the result of a singular eruptive or effusive center from which cryoflows spread out across the more than 400,000 km 2 of Vulcan Planitia. Similarly, the melting-in-place hypothesis is difficult to validate due to energetic considerations and does not account for all of the observations. Instead, we hypothesize that the resurfacing was the result of ammonia-rich cryo-material from the last stages of ocean freezing either buoyantly rising and flowing out on to the pre-Vulcan lowlands, or as a result of more severe disruption that resulted in lithospheric blocks foundering (as suggested by Schenk et al., 2018), and the buoyant, viscous cryo-material under those blocks rising up and spreading out. Under these hypotheses, there would be no singular effusive center, but the sources of the plains material would be in many places across the region, and as the material flowed across the pre-Vulcan lowlands or enveloped the foundering blocks, it would create an extensive plains unit ( Figure 18).
These observations, and our analysis of them, support the conclusion that the extensive plains on Charon are a vast cryoflow emplaced unit (Stern et al., 2015;Moore et al., 2016;, similar to those seen on Ariel and Miranda (Jankowski and Squyres, 1988;Schenk, 1991) and possibly other icy worlds in the solar system. Charon fits into the panoply of icy satellites which display evidence for the movement of cryoflows and resurfacing. ,   to organize and measure the rilles. The Vulcan Planitia region with mapped rilles was divided into 14 areas. These areas were created as individual polygons in a feature class to define each area. The areas overlap their neighboring areas by up to 74 km, and the sizes of these areas (approximately 260 km in diameter) were chosen so that enough rilles were covered by each area to provide meaningful information in the resulting rose diagrams and for statistical analyses. Surface features with topography, like these rilles, are most easily recognized in images with low illumination angles due to the presence of shadows (Binder and McCarthy, 1972;Öhman et al., 2006). The lighting geometry in LORRI images is consistently from the north, which might cause rilles oriented in the E-W directions to be more easily identified. However, it has been found that orientation measurements of topographic features (i.e., impact crater rims) are statistically similar to those taken of the same features with different lighting geometries when the features are manually traced and quantitatively analyzed (Öhman et al., 2006). Since that is our approach here, we find it unlikely that lighting geometry has significantly biased our measured rille orientations.

A C C E P T E D M A N U S C R I P T
The QGIS Line Direction Histogram application was utilized to create rose diagrams that represent the rille azimuth distribution in each of the 14 areas ( Figure 17). To create an accurate azimuth distribution, each recorded azimuth must represent an equal length of the rille. Therefore, to obtain these equallength rille segments, each rille trace was normalized. For each rille segment, the angle and length were calculated, and the angle was used to place that datum in a rose diagram direction bin. The associated length was added to the accumulated length for that bin. The rose diagrams were specified to have 36 direction bins, so that enough data are present within each bin for meaningful analyses. These bins are 5 degrees wide, and centered on whole degrees, so the bin centered on zero degrees covers N-S rilles that have an azimuth from 357.5 • to 2.5 • (and since orientation is really a half-space also covering directions from 177.5 • and 182.5 • , centered on 180 • ).
A C C E P T E D M A N U S C R I P T   The extensive plains on Charon are a result of global expansion and an early subsurface ocean that flowed as a cryoflow and resurfaced a large area.
 The cryoflow most likely was composed of ammonia-rich water ice.
 The cryoflow was similar to Lunar maria emplacement or via a mechanism where lithospheric blocks were severely disrupted.