Tectonism and Enhanced Cryovolcanic Potential Around a Loaded Sputnik Planitia Basin, Pluto

Sputnik Planitia on Pluto is a vast plain consisting of a nitrogen ice deposit filling a broad topographic depression, likely an impact basin. The basin displays a broad, raised rim and is surrounded by numerous extensional fault systems, each with characteristic orientations with respect to the basin center. The nitrogen ice exerts a large mechanical load on the water ice outer shell crust (here also containing the lithosphere). We calculate models of stress and deformation related to this load, varying dimensional, mechanical, and boundary condition properties of the load and Pluto's lithosphere, in order to constrain the conditions that led to the formation of the observed tectonic and topographic signals. We demonstrate that the tectonic configuration is diagnostic of a particular set of conditions that hold for the Sputnik basin and Pluto, including moderate elastic lithosphere thickness (40–75 km, with higher values favored if initial basin topography is compensated) and a basin that was pan‐shaped and shallow (∼3 km) at the time of nitrogen deposition initiation. These tectonic systems show the contributions of both flexural (bending) and membrane (stretching) responses of the lithosphere, with the latter dominating in proportion to the importance of spherical geometry effects. Rim topography may also show an influence of primordial annular trans‐basin ice shell thickening from the impact process. Analysis of stress‐driven cryomagma transport shows that loading stresses can facilitate ascent of cryomagmas in annular zones around the basin, the locations of which overlap the observed distances from Sputnik of several candidate cryovolcanic sites.

• The Sputnik Planitia impact basin on Pluto is filled with nitrogen ice, producing stresses that create outwardradiating fault systems • For likely load distributions, the pattern of faulting is strongly diagnostic of an elastic lithosphere (ice shell) thickness around 50 km • The initial basin depth must not have exceeded several km in order to be consistent with the observed topographic level of nitrogen ice Pluto's encounter hemisphere is dominated by the sprawling Sputnik Planitia (hereafter abbreviated as "SP"), which forms the western portion of the heartshaped, high albedo region of Tombaugh Regio. SP is a massive deposit of predominantly nitrogen ice  measuring ∼1,500 km by ∼900 km (red outline in Figure 1) that partially fills a ∼1,400 km by ∼1,200 km wide elliptical depression termed the Sputnik basin Schenk et al., 2018). Schenk et al. (2018) describe how SP is enclosed by an eroded and modified, broad ridge 125-175 km wide that rises locally up to ∼1 km above surrounding plains and 2.5-3.5 km above the surface of the SP ice sheet, and defined the rim of the basin to be the crest of this broad topographic swell (cyan ellipse in Figure 1); we adopt this definition of the term "basin rim" in this paper. This elliptical basin has been interpreted to form via an oblique impact into Pluto's surface (Johnson et al., 2016;McKinnon et al., 2016;Moore et al., 2016;Nimmo et al., 2016;Schenk et al., 2018), and its initial depth has been estimated to be no deeper than ∼10 km based on gravity scaling the depths of basins on Iapetus . McKinnon et al. (2017) considered the 1 km excess elevation of the basin rim compared to Pluto overall to be consistent with an ejecta blanket. although we note that lateral movement of material due to the collapse of the transient crater also contributes to the structure of the crust around a basin (e.g., Potter et al., 2013), and in some cases the resulting combination of outward flowing crust and ejecta may produce an essentially flat surface rather than exhibiting a positive ("rim-like") topographic expression (see Figures 5-7 of Potter et al., 2013).
The SP basin is thought to be one of the oldest geologic features in Pluto's encounter hemisphere (≳4 Gyr) Schenk et al., 2018), and modeling of volatile behavior in response to topography has shown that infilling of the basin with the majority of surface nitrogen ice would be complete by tens of millions of years after its formation (Bertrand & Forget, 2016;Bertrand et al., 2018;Hamilton et al., 2016), suggesting that the Planitia has been a feature of Pluto's surface for much of its history. The low-viscosity nitrogen ice can be mobilized easily , and where the deposits are interpreted to be thickest toward the center of the Planitia, the cellular morphology of the plains indicate that they are undergoing solid-state convection, powered by the radiogenic heat flow emanating from Pluto's interior Moore et al., 2016;Stern et al., 2015;Trowbridge et al., 2016). The maximum thickness of the deposits filling the basin is not known for certain, but McKinnon et al. (2016) determined that convection cell diameters of 20-40 km imply depths to the base of the nitrogen ice layer of ∼3-6 km, assuming that convection is taking place in the "sluggish lid" regime. This range agrees well with the result of Mills and Montési (2019) that a minimum nitrogen ice load thickness of 4.6 km is required to explain the current topographic profile of the Sputnik basin rim as a flexural bulge forming in response to a thin elastic ice shell being loaded by the nitrogen ice. As indicated by the dashed portion of the cyan ellipse in Figure 1, the SSE portion of the basin rim is absent, with the tapering "tail" of the pear-shaped SP representing overflow of the nitrogen ice from the basin to cover low-lying terrain to the south-indeed, unlike the plains within the central portion of SP, these southern plains do not display a cellular morphology, which has been interpreted to reflect thinning of the nitrogen ice layer to a degree such that it cannot support convection . While erosion may have played a role in explaining the absent southern portion of the rim, it may not have had appreciable relief to begin with owing to the oblique nature of the impact that created the elliptical Sputnik basin, which can result in the omission of rim uplift along the incidence direction for certain impact geometries (e.g., Elbeshausen et al., 2013;Schultz & Gault, 1990).
Reorientation of SP arising from tidal and rotational torques (Keane et al., 2016;Nimmo & Matsuyama, 2007;Rubincam, 2003) can explain the basin's present-day location, but requires the feature to be a positive gravity anomaly, despite its negative topography. The positive mass anomaly associated with the deep basin has been interpreted to be a consequence of the emplacement of several kilometers of dense nitrogen ice on a less dense cooled rigid water ice shell (1.0 and 0.917 g cm −3 , respectively), which is similar to the emplacement of dense mare basalts for some lunar mascon basins (Freed et al., 2014;Melosh et al., 2013), as well as Pluto having a liquid water ocean, with the newly formed basin being initially isostatically compensated by an uplift in the subsurface ocean (which also has a density of 1.0 g cm −3 ) that caused a substantial thinning of the ice shell under the basin (Johnson et al., 2016;Keane et al., 2016;Nimmo et al., 2016). Loading of the basin with nitrogen ice would result in downward deflection of the water ice shell , and Hamilton et al. (2016) offered an alternative hypothesis to explain the basin's formation and location whereby a runaway albedo effect concentrated Pluto's nitrogen ice deposits into a single cap many kilometers thick centered near 30°N (without a pre-existing impact basin being necessary to focus them), with the resulting positive gravity anomaly subsequently locking Sputnik to a longitude directly opposite Charon. Hamilton et al. (2016) argued that the massive accumulation of nitrogen ice onto an early, thin, rigid lithosphere would have caused sufficient downward deflection of the shell to create its own basin. The New Horizons flyby of the Pluto system did not provide any spatially resolved gravitational data and so there is no direct constraint on the compensation state of the Sputnik basin. Several groups have argued that the SP basin has been compensated to some extent (Johnson et al., 2016;Keane et al., 2016;Nimmo et al., 2016). Regardless of how or whether the compensation state of the Sputnik basin has changed, the loading of the basin with this kilometers-thick nitrogen ice deposit has been highly influential for Pluto's subsequent geological history, particularly in terms of governing the configuration of its tectonism.
Pluto's encounter hemisphere displays extensive tectonic deformation in the form of a non-random system of extensional faults (Keane et al., 2016), which indicates global expansion due to partial freezing of a subsurface ocean as the overarching driver of tectonism (Hammond et al., 2016;Nimmo et al., 2016). The variety of configurations and preservation states of the various fault systems suggest multiple deformation episodes and prolonged tectonic activity . The Sputnik basin is centered at 173°E, 25°N on the anti-Charon hemisphere, and close to the Pluto-Charon tidal axis. Keane et al. (2016) found that loading of volatile ices within a Sputnik-sized basin can substantially alter Pluto's inertia tensor, resulting in reorientation of Pluto of around 60° with respect to the rotational and tidal axes, and considered the present location of SP to be the natural consequence of the sequestration of volatile ices within the basin and the resulting reorientation (true polar wander) of Pluto. Finding that tectonism proximal to SP is oriented broadly quasi-radial to SP and that farther away (close to the edge of the encounter hemisphere) is oriented broadly azimuthally, Keane et al. (2016) argued that this orientational transition marks a change in the dominant source of stress, with loading stresses dominating near SP and reorientation stresses dominating farther away. For kilometers-thick deposits within SP, the loading stresses dominate reorientation stresses globally, and the change in stress field with time may be recorded in the crosscutting relationships of the faults, with Keane et al. (2016) predicting that the quasi-azimuthal faults far from SP may be crosscut by the quasi-radial ones closer to it.
Loading of the Sputnik basin may also have created stress conditions in the crust that are favorable for eruption of cryovolcanic material at select locations surrounding it. Tentative cryovolcanic features have been identified on Pluto (Figure 1), including Wright and Piccard Montes, very wide and tall mounds with enormous central depressions at the southern end of SP Schenk et al., 2018;Singer et al., 2016); Hekla Cavus, an oblong, flat-floored depression located within dark uplands west of SP (Ahrens & Chevrier, 2021); a localized zone of mantling identified in and around a segment of Virgil Fossae, interpreted to be ammoniated cryoclastic materials erupted by fountaining events from this segment of the fossae (Cruikshank, Materese et al., 2019;Cruikshank, Umurhan et al., 2019;Dalle Ore et al., 2019); a few fossae and an adjacent impact crater in Viking Terra to the north of Virgil that appear to be infilled with ammoniated, dark material, interpreted to be a water-based cryolava infused with this material that debouched along fault lines and flooded the fossae and crater (Cruikshank et al., 2021); and deep, irregular, flat-floored depressions within Pioneer Terra to the northeast of SP that Howard,  hypothesized may be the source of gaseous emissions from the subsurface, with a methane-rich component depositing on the surrounding terrain to form smooth-textured uplands with broadly rounded divides. Cruikshank, Umurhan et al. (2019) noted that the close spatial association of these various putative cryovolcanic features and SP might indicate a relationship based on enhancement of cryomagma ascent potential in an annular region beyond a large load (i.e., the nitrogen ice filling the Sputnik basin) on Pluto's water ice shell.
Here, we create detailed finite element method (FEM) models of impact-driven lithospheric loading on Pluto (using the COMSOL Multiphysics software package) and evaluate scenarios that are consistent with the spatial distribution and configuration of observed tectonism and proposed cryovolcanic centers.

Structural Mapping and Analysis
Pluto's tectonics have been surveyed previously, less than a year after the flyby (Keane et al., 2016), which was prior to the availability of the latest high quality global mosaic and DTM . The data sets used for our mapping include a global mosaic that covers the encounter hemisphere at a pixel scale ranging from 234 to 835 m/pixel, and the DTM of the encounter hemisphere, which can resolve topographic features as small as ∼1.5 km across and has vertical precision ranging between 90 and 1,120 m . The flyby nature of the New Horizons mission meant that, for the encounter hemisphere, each point on the surface was only imaged at a single solar incidence and emission angle. Assessing topographic relief based on shading in imaging is therefore more difficult in areas around the subsolar point of 130.5°E, 51.5°N, and so the DTM plays a particularly important role in identification of faults here, as well as in areas imaged at oblique angles near the edge of the encounter hemisphere. Figure 2 shows our mapping of tectonic lineations across the encounter hemisphere. We have categorized these lineations into seven distinct classes, each of which appears to bear an orientational and/or stratigraphic relationship to SP. We describe these classes below.

Ridge-Trough System
First described by Schenk et al. (2018), this is a complex, eroded, fragmentary, NNE-SSW-trending band of graben, troughs, ridges, plateaus, tilted blocks, and elongate depressions (the latter occurring within Cthulhu Macula and mapped as RTS Depressions in Figure 2) that measures ∼300-400 km wide and extends at least 3,200 km from the north pole to the limit of coverage at ∼45°S, crossing the equator at ∼155°E. The structure may well extend further into the poorly resolved far side, and into the shadowed southern regions Stern et al., 2021). The system certainly represents the earliest evidence of tectonism yet seen on Pluto due to its highly eroded state, and because its elements are invariably crosscut by other tectonic lineations. It may even predate the Sputnik basin-forming impact, but since it crosscuts the broad raised rim of the Sputnik basin and terrain leading down to SP itself, some deformation did still occur after the basin formed . Equatorial crustal thickening has been hypothesized to be the cause of such an immense tectonic feature aligned along a great circle, although this would require the system to be aligned along a "paleo-equator" prior to reorientation of Pluto (McGovern et al., 2019).

Sputnik-Azimuthal
In the far west and far east of the encounter hemisphere lie Djanggawul and Mwindo Fossae respectively, tectonic systems that are oriented circumferentially about a pole located at 170°E, 25°N, very close to the geographic center of SP, and which implies an origin tied to that feature. The faults are removed from the pole by ∼1,500 and ∼1,100 km, respectively. Mwindo Fossae are unusual in that they converge to a nexus, implying that a localized stress field caused them to diverge from the Sputnik-azimuthal orientation (McGovern et al., 2019). The geology at the nexus is unremarkable relative to its surroundings, providing no clue as to what endogenic process may have contributed to this convergence. The longest fault of Mwindo Fossae, named Sleipnir Fossa, extends 560 km to the SW from the nexus. The modeling of Keane et al. (2016) found that reorientation of Pluto in response to infilling of the Sputnik basin with nitrogen ice would generate stresses that are approximately consistent with the Sputnik-azimuthal orientations of these systems, and that the azimuthal faults might be crosscut by the quasi-radial ones proximal to SP (mapped as the East and West Sputnik systems in Figure 2). We find that these azimuthal and radial systems are sufficiently removed from each other, however, such that they do not intersect, meaning that crosscutting relationships cannot be established. But the azimuthal faults do appear older than the radial ones, in that they are more eroded and are superposed by impact craters in a few places (the radial faults, particularly those of West Sputnik, always crosscut craters that they encounter). We note, however, that the less well-preserved appearance of the azimuthal faults may in part be a consequence of them being located near the edge of the encounter hemisphere, where the pixel scale of New Horizons imaging is coarse (∼475 m/pixel) and emission 10.1029/2021JE006964 5 of 31 angles are high (>70°), which can make the faults appear less sharp relative to those nearer to SP and the center of the encounter hemisphere (covered by imaging with pixel scale of ∼315 m/pixel and emission angles of 30°-60°).

West Sputnik
The terrain to the west of SP (between 100° and 150°E) displays the best-preserved tectonism in Pluto's encounter hemisphere. The West Sputnik system consists of generally sharply defined graben and troughs, which at their northern extent are quasi-radial to SP (Inanna and Dumuzi Fossae are prominent examples of these), but which bend southwards to adopt a N-S orientation in the southern hemisphere (Hermod Fossae). Between 3°N and 5°S, these faults intersect with the NE-SW-aligned faults of the Sputnik-Radial system (Virgil Fossae and Beatrice Fossa). The West Sputnik faults crosscut all craters that they encounter, indicating their relative youth. The modeling of Keane et al. (2016) indicated that these and other quasi-radial fault systems to the east and west of SP (Sputnik-Radial, Northwest Sputnik and East Sputnik as mapped in our study) might have formed in response to loading of the Sputnik basin with nitrogen ice.

Sputnik-Radial
This fault system consists of sharply defined graben and troughs that are similarly well preserved as those of West Sputnik, and which crosscut all craters they encounter. The main components of this system are Virgil and Beatrice Fossae which, along with a similarly oriented fault located at 30°S 120°E, form a belt with a consistent NE-SW orientation that is essentially radial to the center of SP. It is the consistently Sputnik-radial orientation of these faults across their entire observable lengths that differentiates them from those of the West Sputnik system, which display a quasi-Sputnik-radial orientation proximal to SP, but start to bend southwards a few hundred km from the edge of SP. At 8°N, 126°E (indicated by the white asterisk in Figure 2a), components of Virgil Fossae form what has been interpreted to be a set of strike-slip duplexes caused by dilatational dip-slip normal faults with a right-lateral strike-slip component (Cruikshank, Umurhan et al., 2019). Faults belonging to Virgil and Beatrice Fossae appear to crosscut those of Hermod Fossae and vice-versa, suggesting that, despite their different orientations, the faults of the West Sputnik and Sputnik-Radial systems likely formed concurrently. Ammoniated water ice deposits that appear to mantle underlying terrain have been identified in and around a segment of Virgil Fossae (Dalle Ore et al., 2019), possibly representing cryoclastic materials recently (<1 Gyr) erupted by fountaining events from this segment of the fossae (Cruikshank, Materese et al., 2019;Cruikshank, Umurhan et al., 2019).

Northwest Sputnik
At its southern extent, the faults of this system appear to be almost continuous with those of the West Sputnik and Sputnik-Radial systems. But as the Northwest Sputnik system is traced farther northwards, the divergence of the orientation of its faults from those of West Sputnik becomes increasingly marked, such that at the northern extents of both systems their faults are oriented at an oblique angle to one other. The convergence of all three of these quasi-Sputnik-radial systems at ∼20°N 142°E may indicate that they formed during the same tectonic episode brought on by lithospheric stresses related to the formation of the Sputnik basin and its infilling with nitrogen ice (Keane et al., 2016). The Northwest Sputnik faults are quasi-radial to SP at the southern extent of the system, becoming radial at the northern extent. They tend to be smaller scale than those of the West Sputnik and Sputnik-Radial systems, generally forming short, narrow troughs and scarps rather than wide graben 100°km long. Where they intersect with the West Sputnik faults, it is quite difficult to determine crosscutting relationships due to the narrowness of the Northwest Sputnik faults: in some cases those of Northwest Sputnik seem to crosscut those of West Sputnik and vice-versa, which may be regarded as evidence for the roughly contemporaneous formation of the two systems. Cruikshank et al. (2021) identified morphological evidence of infilling of the Uncama Fossa graben (in the southern part of the Northwest Sputnik system), as well as an adjacent impact crater, with ammoniated, dark material. They suggested that the crater and fossa trough might have been recently (∼1 Gyr) flooded by a cryolava debouched along fault lines in the trough and in the floor of the impact crater.

East Sputnik
Tectonism in the area to the immediate east of SP (between 185°E and 225°E) exhibits a diverse range of morphologies, including pit chains, narrow troughs, tall scarps with scalloped walls, and broad graben reaching a few tens of km across. The alignment of the pit chains parallel to structural trends in the vicinity suggests that they are where surface collapse has occurred as tectonism disturbs an overlying mantle . Indeed, this region incorporates a number of terrains that are interpreted to have experienced largescale deposition of methane-rich material since the formation of the Sputnik basin, including the bladed terrain deposits to the east  (and also the bright, pitted uplands separating SP from the bladed terrain, interpreted to be a modified, westerly extension of the bladed terrain deposits) and the smooth uplands to the northeast . The faults of the East Sputnik system are predominantly oriented NW-SE, with those in the south being quasi-radial to SP, whereas those in the north are oriented obliquely to radial trends, and even to nearly tangential orientations near the northeast basin margin. Great circles extrapolated from the northern faults of the East Sputnik system align fairly well with the northern faults of the West Sputnik system, raising the possibility that these geographically separate systems may have formed as a consequence of a single tectonic episode, with lithospheric stresses being mirrored on either side of SP. The East Sputnik faults generally have a more degraded appearance that those of West Sputnik and in a few instances are superposed by impact craters. Their degraded appearance is at least partly due to the major mantling and erosional episodes that have affected this region, which the faults west of SP have not been subjected to. In addition, the fact that the faults of East Sputnik often manifest as pit chains (not seen elsewhere in the encounter hemisphere) and are currently experiencing ongoing deposition of volatile ices where they occur in the bright, pitted uplands, further complicates assessment of their relative age based on preservation state.

Southwest Sputnik
This is the smallest and most tenuously defined tectonic system, which features localized clusters of NE-SWtrending ridges, scarps, and troughs on the southwestern rim of the Sputnik basin. Their sparseness lends uncertainty to the interpretation that they share an origin in a single episode of tectonism, although their quasi-radial orientation to the center of SP raises the possibility that they originated through lithospheric stresses associated with the basin.

Undesignated
All encounter hemisphere tectonics that are not classified within one of the seven aforementioned systems are termed "undesignated" for the purpose of this study. The orientation of these faults bears no obvious stratigraphic or (quasi-)radial/azimuthal relation to the center of SP, and some may originate from localized crustal stress conditions, for example, a cluster of reticulate networks of fractures that occur to the west of Wright and Piccard Montes, and the series of E-W-oriented faults that form the non-Sputnik-azimuthal elements of Mwindo Fossae. Other systems, however, are regional in their extent, including a series of E-W-and WNW-ESE-trending faults that appear to crosscut Djanggawul Fossae, and a sparsely populated belt of WNW-ESE-trending faults (including Kaknú Fossa) that crosses Hermod Fossae at 30°S 120°E.

Topography of Basin-Filling Units and Surroundings
A global stereo digital elevation model (DEM) for Pluto was created from images collected by the MVIC and LORRI instruments . The topographic character of both the interior of the Sputnik basin and of the terrain surrounding it is a key input to creating mechanical models of SP's loading of Pluto's icy shell lithosphere. To facilitate model generation, we take profiles across this DEM that originate from a nominal center of symmetry of SP (located at 25°N 173°E) for sectors representative of the axisymmetric nature of the northern part of SP. Profile azimuths were selected to avoid strongly non-axisymmetric influences such as the south-southeast extending embayment of SP and the ridge-trough system (RTS; purple lines in Figures 1 and 3). Individual profiles within the eastern group ( Figure 4a) show maximum relief from rim to lowest point (at the outer margin of the SP nitrogen ice deposits) approaching almost 4 km. However, the mean profile shows 3 km for this difference. Additionally, the topography surrounding the basin is on average lower than the rim crest by about 0.8 km. We have subtracted a reference value of the far-field (radius r > 700 km) topography, 0.7 km elevation, from the mean profile to arrive at an average elevation of the interior nitrogen ice plains of −2 km; the plains tend to be slightly higher than this value near the basin center and lower near the margin.
The western group of profiles ( Figure 4b) shows lesser topographic variation exterior to the basin rim than those of the eastern group ( Figure 4a). We subtract the same reference far-field topography value (0.7 km) from the former as from the latter, and find an ice-margin-to-rim-peak difference approaching 3.7 km (Figure 4b). The western mean referenced elevation profile again shows an elevation slightly higher than −2 km near the center and closer to −3 km at the margin. This transition appears to occur at a much higher value of r than for the eastern mean referenced profile, due to anomalously high elevations in the r = 300-400 km range. However, individual profiles that avoid the al-Idrisi, Zheng-He, and Baret Montes, (which are attributed to water ice blocks calving off the basin margins; White et al., 2017;O'Hara and Dombard, 2021;Skjetne et al., 2021) show a mid-basin dip below −2 km, consistent with the eastern profile. Since we consider the blocks in the various Montes to be "anomalous" and therefore not representative of the character of the general SP nitrogen ice load, we adopt a −2 km offset between the basin surface and the far-field ice shell surface in the models presented below.

Modeling Method
In order to constrain the mechanisms by which prominent tectonic systems surrounding SP have formed, we seek to explore the range of possible lithospheric responses (stress and deformation) to loads infilling the Sputnik basin. To this end, we use the COMSOL Multiphysics FEM code (Version 5.6) to calculate models of the elastic response of Pluto's icy shell lithosphere to infill of a Sputnik-sized impact basin (radius ≈ 500 km) by nitrogen ice. For computational economy, we use 2-D spherical axisymmetric geometry, with a shell encompassing the entire circumference of the planet. The nominal elastic shell (lithosphere) thickness T e is set to 50 km, with larger and smaller values tested as well. Local shell thickness variations are created corresponding to the initial topography of the Sputnik Planitia impact basin, and in some cases, compensation of the basin topography applied as relief on the ocean-shell boundary (or "OSB"); in such cases, T e represents the far-field constant value of shell thickness. We utilize the "free triangular" meshing scheme of COMSOL Multiphysics to create a finite element mesh within stated elastic shell boundaries (Figures S1 and S2 in Supporting Information S1). The radial distance    for Sputnik Planitia (SP; colored dashed lines), drawn from the geometric center of the SP basin, taken to be 25° N 175° E, sampled at 5 km increments. Horizontal axis: along-surface distance r from SP center. Left vertical axis: elevation above mean planetary datum; right vertical axis: local slope (degrees) calculated from adjacent points on profile. The black solid line denotes the average of the profiles. (a) Eastern profiles, at azimuths of 45°, 60°, 75°, 90°, and 105°. (b) Western profiles, at azimuths of 240°, 255°, 270°, 285°, and 300°. The average elevation of the terrain surrounding the basin, taken as 0.7 km for both profile groups, was subtracted from the raw topographic profiles to give relative elevation, which is the most suitable basis for comparing to the finite element method models. The locations of prominent topographic features are indicated in red. Note that the vertical exaggeration of the topography in both panels (a) and (b) is 184×. The long-baseline slope from the edge of SP to the basin rim crest ranges between 1° and 2°. of the top surface of the shell from the origin corresponds to the radius of Pluto (R p = 1,188 km), adjusted locally to account for the preexisting basin topography (see below). The bottom surface corresponds to a radius R p -T e . Individual elements are constrained in size by a maximum limit T e /20, yielding elements of order 1-2 km in size. We plot model quantities using distance reckoned on the surface of Pluto at R p = 1,188 km, termed the "projected radial distance" r proj as described in the Supporting Information S1.

Basin Topography and Load Dimensions
The magnitude and shape of the infill load is constrained by topographic data for the central basin and surrounding uplands ( Figure S1 in Supporting Information S1) and by insights from hydrocode models of impacts into planetary bodies such as Pluto (Johnson et al., 2016) and the Moon (Freed et al., 2014;Melosh et al., 2013;Potter et al., 2012Potter et al., , 2013. Ultimately, it will reflect to some extent the shape of the basin at the time of initiation of infill. To represent the initial basin shape, we use a "super-Gaussian" profile of the type used to model the topography of oceanic lithospheric swells on Earth (Wessel, 1993), where r proj is projected radial distance as defined above, d bc is the depth at the center of the basin (on the model symmetry axis, r proj = 0), w c is a characteristic width, and p is an exponent equal to 2 for a standard Gaussian profile (in which case w c is the half-width), with increasing values of p producing a zone of flat topography of increasing width emanating from the basin center. This negative relief profile is applied to the upper surface of the model shell as topography. We assign w c to give a basin that reaches a given r proj at 1% of the central depth of the basin. As a baseline, we use a flat-bottomed basin profile with super-Gaussian exponent p = 6 to resemble observations of relatively "fresh" or "pristine" basins (e.g., Potter et al., 2013) and the results of hydrocode impact models (e.g., Johnson et al., 2016;Potter et al., 2012). We also test "bowl-shaped" basins with a standard Gaussian (p = 2) profile and "hyper-flat" (essentially disk-shaped) basins with p = 10. Given the initial basin shape, we can calculate the magnitude of the initial stress state in the shell, here taken to be lithostatic: each normal stress component is set equal to ρ c g (z-H b ).
We apply nitrogen ice basin-filling loads as surface force boundary conditions, using an iterative process to determine load configurations that are consistent with observations. Successful models will produce ∼2 km relief between the load surface and the level of the surrounding plains ( Figure 4). This topography-matching requirement, in combination with the initial basin shape and shell deflection profile, will determine the ultimate thickness and shape of the load, thereby providing a constraint on nitrogen ice volumes within the basin.
The basal boundary condition at the bottom of lithosphere comprises two parts: the first is a restoring pressure proportional to the basal deflection d (sometimes called a "Winkler" foundation) of magnitude ρ o g d, where ρ o is the density of ocean, reflecting the buoyant support from the ocean; the second is a "counterweight" to the initial state of lithostatic stress enforced within the lithosphere (e.g., Galgana et al., 2013Galgana et al., , 2011Le Corvec et al., 2015), of magnitude ρ c g t c , where t c is the local thickness of the shell (accounting for the initial basin topography and any compensation), establishing equilibrium of the unloaded shell.
Hydrocode models of planetary basin-forming impacts also predict significant topography at the basal crustal boundary (e.g., Johnson et al., 2016;Melosh et al., 2013;Potter et al., 2013Potter et al., , 2012. This topography creates post-impact lithospheric uplift that plays a strong role in creating the gravity signature of lunar mascons (e.g., Andrews-Hanna, 2013;Freed et al., 2014;Melosh et al., 2013). To represent such effects, our models also include a "crustal collar" buoyant load at the base of the lithosphere, characterized with a simple Gaussian shape. For the nominal model we use a Gaussian half-width w cc = 65 km and a load center at r plc = 600 km, reflecting crustal thickening expected from the impact process. Both w cc and r plc are varied to test different relationships of the crustal collar to basin topography. We note that when we refer to a "crustal collar" load we are referring to the subsurface structure and not any topography that the subsequent response generates, and also not to any surface topography created by ejecta.

Basin Compensation State
The physical structure of Pluto's ice shell is interlinked with the compensation state of the basin topography, which itself is bound to the thermal state of the shell and basal ocean. In one class of thermal models of Pluto, the lowermost parts of the shell are thought to be warm enough to undergo ductile deformation through creep processes in the ice (e.g., Bierson et al., 2020Bierson et al., , 2018Kimura & Kamata, 2020). A second class of Pluto thermal models posits a uniformly cold shell with no ductile lower layer, likely enabled by a layer of insulating material such as clathrate hydrates (Kamata et al., 2019). The former class of models would tend to erase compensating OSB relief over time via ductile flow, while the latter could preserve it until the present day. However, the former type of models are unstable with respect to surface loading. Flexure-induced OSB relief supports the load through a buoyancy effect: if this relief is removed the plate would both lose support and get mechanically thinner, resulting in more deflection, which would in turn generate more relief, followed by more removal of relief, more deflection, and so on. In other words, this scenario generates runaway deflection, an instability that would render support of a surface feature impossible. This means that preservation of the OSB relief is an inherent assumption of our models, and that therefore we are implicitly assuming a cold insulated shell scenario such as that proposed by Kamata et al. (2019).
We take two main approaches with respect to compensation state, constituting end-member scenarios intended to bracket the actual structure of Pluto's icy shell at SP. In the first, we do not apply any initial relief at the basal OSB interface. Such a case corresponds (excepting basin relief) to icy shell loading model solutions that assume a constant-thickness shell (e.g., Janes and Melosh, 1990;Keane et al., 2016). In the second approach, compensation is accounted for by relief on the OSB via an Airy-type model: . For such cases, T e is the far-field shell thickness, but the sub-basin shell will be locally thinned . Note that the compensation magnitude (the term in parentheses) equals 11.5 for Pluto, meaning that the shell thickness at the basin center will vanish for sufficiently low values of T e and/or high values of d bc . In general, the immediate post-impact mechanical state of an impact basin is not perfect isostatic equilibrium, and post-impact adjustments based on thermal and mechanical responses can result in subsequent changes in surface topography and internal layer relief (e.g., Johnson et al., 2016;Melosh et al., 2013;Potter et al., 2013Potter et al., , 2012. Such responses might include ductile flow of the impact-warmed sub-basin portion of the shell, potentially removing some or all of the initial relief on the surface and OSB. For this reason, we consider our compensated and uncompensated models to represent end-members of the potential state of relief on the OSB. Finally, the "crustal collar" loads described above address one of the principal departures from simple Airy-type isostatic compensation seen in impact basin settings. Given that the response time to post-impact dis-equilibrium is set by the viscosity of the sub-crustal material (in this case, Pluto's ocean), the response to a crusal collar load would be essentially instantaneous, thereby setting the initial stress state for the region.

Faulting Regimes
We characterize the faulting type predicted by the stress tensor within the shell using the Aψ parameter (Simpson, 1997). Values range over ±180°, with specific fault types corresponding to the labels above the brightest colors in Figures 5-12 at values ± 150° (thrust), ±90° (strike-slip), and ±30° (normal), with the sign determining the specific orientations of the faults, as labeled in the figures. Values in-between these represent mixed modes of faulting. Values of ±180° correspond to compression with ambiguous orientation (termed "pure constriction"), and similarly for value of 0° for extension ("pure extension"). We will identify the presence of a specific fault type plus orientation prediction as a stress "regime" (e.g., radial normal regime zone) rather than say "radial normal fault zone," because faulting per se is not predicted to happen where the failure criterion is not satisfied.

Material Properties and Failure Criterion
We adopt material properties appropriate to the water ice and nitrogen ice constituents of the models (Table 1). We recognize that strength-type properties of deformed large-scale rock (ice) assemblages often fall short of laboratory derived values. Following this philosophy, we adopt a Young's Modulus value for water ice of 5 MPa, intermediate between values of order 1 MPa from field observations (e.g., Vaughan, 1995) and 9 MPa from laboratory-scale specimens (e.g., Petrenko & Whitworth, 1999); See also the discussion in Nimmo (2004). We construct Mohr-Coulomb failure envelopes with parameters cohesion c = 1 MPa and angle of internal friction angle ϕ = 30°.

Stress Definitions
We adopt an axisymmetric spherical shell system with three primary normal stress components. The locally horizontal components are σ h and σ ϕ , in the plane of the model and perpendicular to it, respectively; σ ϕ is sometimes called the "hoop stress." The vertical component (i.e., pointing to the central point contained within the shell) is called σ v . We also define a "tectonic stress" after Rubin (1995) for each horizontal component: These components are useful to characterize the tectonic implications of specific stress tensor configurations and for calculating cryomagma ascent criteria within dikes. We use the "engineering" convention that tension is positive and compression is negative in sign.

Cryomagma Ascent Criteria
Following McGovern et al. (2013McGovern et al. ( , 2016, we use two criteria for cryomagma ascent: (a) the stress orientation criterion requires that the least compressive stress be oriented horizontally to allow vertical dikes to form (e.g., Anderson, 1951). This criterion can be expressed as σ T > 0, for either of the tectonic stress components. (b) We use the formulation of Rubin (1995) to calculate cryomagma ascent velocity u z within vertical dikes, using the vertical gradients of tectonic stress components dΔσ T /dz, where η is (cryo)magma viscosity, w is dike width, Δρ is density contrast between magma and host rock (here, between liquid water and the water ice shell), g is planetary gravity, and ΔP is magma overpressure. The first and second terms in parentheses on the right hand side of Equation 2 can be equated to yield a solution for an effective buoyant density Δρ efb , that can counteract the negative buoyancy of water in ice (here taken to be −80 kg/m 3 ).

Model Nomenclature
To facilitate concise and efficient identification of the various models presented herein, we constructed a list of model labels (Table 2). Each label contains information on the compensation state, T e , the presence or absence of crustal collar loads, and variations in basin depth and shape in a compact 4 or 5 character format.

Modeling Results
We first describe the stress state and deformation ( Figure 5) induced by a baseline ("nominal") model U50N with an initial "pan-shaped" basin with d bc = 3 km, p = 6, and r proj = 500 km (yielding w c = 387.6 km), and without OSB compensation. Near the symmetry axis, the applied basin-filling load produces horizontal stress components σ h and σ ϕ that are compressional throughout almost the entire thickness of the lithosphere. At the surface of Note. Material properties adopted for the finite element method models. Note. Model names consist of a 1-letter prefix denoting the compensation state of the initial basin topography, 2 digits denoting the far-field elastic thickness of the shell T e in km, and a suffix of up to 3 characters that denotes information about basin shape and presence or absence of a crustal collar load. Note that all models with a suffix beginning with "N" (representing the word "Nominal") have d b = 3 km and a pan-shaped initial basin.

Table 2
Model Nomenclature the lithosphere, the magnitudes of these stresses (red and black lines in Figure 5a, respectively) are nearly identical (with |σ ϕ | slightly higher in mid-basin), resulting in a fault regime prediction of generic compression ("pure constriction," light purple colors in Figure 5), although any such faults would be obscured by the nitrogen ice load. In the outer region of the basin (380 km <r proj <4,300 km), with increasing radial distance r proj , the stress regimes cycle through narrow zones of radial thrust, strike-slip, and concentric normal, ultimately reaching a wider region characterized by the radial normal regime out to 640 km, as σ h and σ ϕ diverge in magnitude. Further outward, a strike-slip regime extends to well beyond r proj = 1,000 km.
The failure criterion is exceeded in two surface regions: a region beneath the load (0<r proj <270 km) characterized by pure constriction, and a broader zone beyond the load (450<r proj <850 km) that comprises proximal radial normal faulting and distal strike-slip faulting (Figure 5b). For r proj >640 km, a strike-slip regime is seen at the bottom and through most of the depth of the lithosphere, and for r proj > 650 km the entire lithosphere is in strikeslip mode ( Figure 5). These findings stem from the extensional out-of-plane stress σ ϕ produced by the membrane response of a curved (spherical) lithosphere (Turcotte et al., 1981) The effective buoyant density calculated from the vertical gradient of the out-of-plane tectonic stress Δρ efb is elevated, with peak values of about 60 kg/m 3 , in a region 450 <r proj <650 km, providing a partial offset of the negative buoyancy of water in ice. This offset appears strongest at the bottom of the lithosphere. This region also satisfies the stress orientation criterion (σ ϕ > 0) throughout its radial width, owing to the extensional membrane contribution to the out-of-plane stress σ ϕ throughout the thickness of the lithosphere (see Section 6.2). . Regions for which σ ϕ is not the most extensional stress (i.e., where the stress orientation criterion for radial diking is not met) are masked out by gray.
Adding a subsurface buoyant load ("crustal collar") reflecting the post-impact configuration of the shell can significantly change the shell stress state, with important implications for stress regimes, observed faulting, and cryomagma ascent potential ( Figure 5). For a model with crustal collar characteristic width 65 km and center projected distance 650 km (model U50CO, Figure 6), the lateral and vertical extents of the radial normal regime zone are larger relative to the model U50N in Figure 5, as are those of the region that satisfies the failure criterion.
Thus, the addition of the crustal collar load significantly enhances the potential for normal faulting oriented radially to Sputnik basin, although we note that the Aψ parameter moves close to pure extension near the center of the zone. The effective buoyant density calculated from the vertical gradient of the out-of-plane tectonic stress shows a region ranging from ≈580 to 720 km where Δρ efb is greater than 80 kg/m 3 (maximum value 139 kg/m 3 ), thereby allowing cryomagma ascent despite the negative buoyancy of water in ice. This region also satisfies the stress orientation criterion (σ ϕ > 0) throughout its entirety, again owing to the extensional membrane component of stress.
Lithospheric thickness plays a key role in modulating the stress response to loads. For a model with shell thickness T e = 30 km (model U30N, Figure S3 in Supporting Information S1), the configuration of stress is markedly different from model U50N, with a large separation of the horizontal normal stress magnitudes at the surface ( Figure S3a in Supporting Information S1) and a prediction of strike-slip regime throughout almost the entire lithosphere beyond the basin. Only a tiny excursion of σ h above the vertical normal stress (= 0 beyond the load) at the surface for r proj between ∼420 and 540 km produces a narrow and extremely shallow zone of radial normal regime. The depths of predicted failure regions beneath the load and beyond it are both substantially deeper and broader than for the nominal case ( Figure 5), nearly penetrating the entire thickness of the shell ( Figure S3b in Supporting Information S1). This condition reflects the greatly elevated stress magnitudes resulting from the greater load magnitude, as required to meet the topographic constraint on the weaker lithosphere, see Figure 7. Cryomagma ascent is greatly enhanced in an annular zone 450 km <r proj <600 km, with Δρ efb values well in excess of 80 kg/m 3 .
For a T e = 30 km case with a crustal collar model (model U30CO, Figure S4 in Supporting Information S1), the stress regime is broadly similar to the model without the collar, but with the diminishment of the small radial     Figure  S3 in Supporting Information S1). Note that white color of outer dashed arc denotes an outer limit of predicted failure that exceeds the bounds of the figure. (b) Model U50C with T e = 50 km ( Figure 5). (C) Model U70C with T e = 70 km ( Figure S5 in Supporting Information S1).
normal regime zone near the basin margin and the creation of a larger radial normal regime zone above the collar load (580 <r proj <710 km), and a band with Aψ values approaching 60° between these two zones. Similarly, the zone of enhanced Δρ efb moves outwards in r proj to reside above the collar load.
A model with T e = 70 km (model U70N, Figure S5 in Supporting Information S1) shows substantial differences from the models with thinner lithospheres, but also some similarities. As with model U50N (Figure 5), when T e = 70 km the horizontal stresses at the surface of the shell are similar in magnitude to each other (although Figure 10. As in Figure 9, for uncompensated models with crustal collar loads. (a) Model U30CO with T e = 30 km ( Figure  S4 in Supporting Information S1). Note white outer dashed failure region boundary. (b) Model U50CO with T e = 50 km ( Figure 6). (c) Model U70CO with T e = 70 km ( Figure S6 in Supporting Information S1). somewhat lower than their equivalents when T e = 50 km). However, the surface stress regimes at the outer part of the basin for model U70N go through a cycle of radial thrust, strike slip, concentric normal, and radial normal between 360 and 730 km, with only the latter two potentially visible beyond the edge of the load (r proj ∼ 450 km). The zones of predicted failure are markedly narrower in r proj and shallower in depth than for model U50N, with the normal fault zone being barely 80 km in extent and entirely concentric in character. Peak Δρ efb values of about 15-20 kg/m 3 are also much smaller than for model U50N. The latter properties of this model reflect the much thinner load compared to the models with thinner shells, again dictated by the topography constraint (Figure 7).
A model with T e = 70 km and a crustal collar load (model U70CO, Figure S6 in Supporting Information S1) shows a shrinking of the surface strike-slip regime zone, and a significant expansion of extent and depth of the zone containing the two normal fault regimes, relative to the results for models U30CO and U50CO. Of these, the concentric normal regime zone covers most of the surface (440 km <r proj <690 km), although the most distal part of this combined zone, and most of its depth, fall into the radial normal regime Values of Δρ efb approach 50 kg/ m 3 above the crustal collar ( Figure S6c in Supporting Information S1), comparable to those observed for model U50N which lacks a collar load (Figure 5c).
To investigate the effects of potential variations in the pre-loading central (maximum) depth of the basin d bc , we calculated pan-shaped basin models with d bc = 4, 5, and 6 km (and otherwise nominal conditions). The results indicate that merely doubling d bc to 6 km (model U50D6) almost quadruples the final maximum post-loading thickness of the nitrogen ice load d final , (29 vs. 7 km for the model U50N d bc = 3 km case: Figure 7b) demonstrating a highly non-linear relationship between d bc and d final . The stress state produced for such a model (U50D5, d bc = 5km, Figure S7 in Supporting Information S1) resembles that for the T e = 30 km case (U30N, Figure S3 in Supporting Information S1) in terms of stress magnitudes ( Figure S7a in Supporting Information S1), fault type regimes, ( Figure S7b in Supporting Information S1) and the extreme depth of the failed region (Figures S7b and S7c in Supporting Information S1), in this case penetrating the entire shell for 450 km <r proj <510 km. Cryomagmatic potential is greatly enhanced, with the critical 80 kg/m 3 value of Δρ efb exceeded throughout the entire lithosphere for 520 km <r proj <710 km.
We also consider models with an initially "bowl-shaped" basin (standard "Gaussian" shape with exponent p = 2). Such a model (U50B3) with otherwise nominal parameter values, including initial basin center depth d bc value of 3 km, exhibits very small stress magnitudes, because the central concentration of the bowl shape results in very small load volume (Figure 7). For a model with d bc = 6 km (U50B6, Figure S8 in Supporting Information S1), the results resemble models U70N and U70CO) in that a significant zone of concentric normal faulting regime is predicted at the surface. However, in contrast to the T e = 70 km models, this zone occurs within a wide zone of failed shell (within the black contour on Figure S8b in Supporting Information S1), thus yielding a prediction of actual concentric normal faults being visible at the surface, a prediction that contradicts observations. The trans-basin zone of enhanced Δρ efb values greater than the critical density offset of 80 kg/m 3 takes on an unusual "funnel" shape, wider at the bottom of the lithosphere (∼140 km, see Figure S8c in Supporting Information S1).
For a given T e , models in which the initial basin topography is compensated exhibit substantial departures from the results of the corresponding uncompensated models, stemming from reductions in sub-basin shell thickness that can be substantial for lower values of T e and complete (the central shell vanishes) for T e <37.5 km with the nominal value of d bc (3 km). For T e = 50 km (model C50N), the shell thickness at the symmetry axis is 12.5 km (or T e /4), greatly reducing the ability of the sub-basin shell to support the basin-filling load. Most of the shell between the symmetry axis and r proj = 360 km experiences stress levels well in excess of the failure criterion (Figure 8), and the peak magnitude of deflection ( Figure S11 in Supporting Information S1) is more than double the amount for the uncompensated case U50N (compare with Figure 7a), creating a load volume 1.65 times that of the nominal model Trans-basin surface faulting is predicted to be almost entirely in the strike-slip regime ( Figure 8 and Figure S10a in Supporting Information S1), in contrast to the dominantly concentric normal regime seen close to the basin for the equivalent nominal model C50N ( Figures 5 and 9b), but similar to the thinner lithosphere model U30N case ( Figure S3 in Supporting Information S1 and Figure 9a). Cryomagmatic potential is greatly enhanced in the upper lithosphere at the margin of the basin (Δρ efb >> 80 kg/m 3 , Figure 8), resembling the magnitudes seen for model U30N km ( Figure S3 in Supporting Information S1 and Figure 9a), but the principal stress orientations in the lower lithosphere are different, causing a discontinuity in the orientations of any dikes that might form.
For a compensated model with T e = 70 km (C70N, Figure S9 in Supporting Information S1), the central shell thickness is 37.5 km (or slightly greater than T e /2). Compared to the equivalent uncompensated model U70N ( Figure S5 in Supporting Information S1 and Figure 9c), the region of predicted failure beyond the basin has greater radial extent while the extent of the concentric normal faulting regime is reduced ( Figure S9 in Supporting Information S1 and Figure 10c), such that it is covered by the load and would not be observed at the surface, and the overall pattern most closely resembles that of model U50N ( Figure S5 in Supporting Information S1 and Figure 9b). The zone of enhanced Δρ efb is again located right at the margin of the basin ( Figure S9 in Supporting Information S1), but with substantially reduced magnitude compared to the equivalent uncompensated model ( Figure S5 in Supporting Information S1 and Figure 9c).
Models with basin compensation and crustal collar loads also show significant differences from the corresponding uncompensated models at the same T e . Such a model with T e = 70 km (C70CO, Figure S10 in Supporting Information S1) exhibits a broad region beyond the basin margin in a radial normal stress regime, with the failure criterion satisfied for 390 km < r proj <841 km. Such a configuration strongly resembles that of the uncompensated crustal collar model at T e = 50 km (U50CO, Figure 6), and is markedly different from the uncompensated crustal collar model at T e = 70 km (U70CO, Figure S6 in Supporting Information S1). The zone of enhanced Δρ efb has similar magnitude and location to that for the equivalent uncompensated model (U70CO, Figure S6 in Supporting Information S1).

Comparison of Model Results to Observed Tectonics
Our models give insights into the lithospheric conditions on Pluto at the time of loading that produce the observed distribution of faulting, fracturing, topography, and cryovolcanism in the vicinity of the Sputnik basin. As stated in Section 2, the fault systems surrounding SP are characterized by a generally extensional nature and predominantly radial orientations near the rim/margins of the topographic basin, with a tendency toward more oblique orientations with increasing distance from the center of SP. To facilitate the evaluation of successful and unsuccessful model scenarios for explaining the observed tectonic systems surrounding SP, we plot the Aψ parameter and extent of predicted failure for several models side-by-side (Figures 9-12). We find that for each class of basin compensation, there is a best-fit T e value that produces stress magnitude and orientation predictions that match the observed fault distributions at SP well, with values significantly above and below those values producing poor fits to the observations. Basin compensation decreases the thickness of the shell beneath the basin, making it respond to basin-filling loads as an overall thinner shell. Thus, we find that the best-fit value of T e is somewhat higher for compensated models than uncompensated ones.
The specific fault systems identified in our mapping (Figures 2 and 3) can be used to evaluate the effectiveness of our various models in representing the conditions that controlled faulting on Pluto. For example, the "Sputnik-Radial" system (green in Figures 1 and 14), extending from r proj ∼ 600-1,600 km, exhibits dominantly radial-to-SP orientations (hence the name) with somewhat sinusoidal variations around a baseline radial trend. Such orientations in the observed locations are most consistent with model surface conditions beyond the basin margin that exhibit a radial normal faulting regime (Aψ parameter) and satisfy the failure criterion. For models lacking compensation of initial basin topography ("U" prefix), such conditions are produced for T e = 50 km, such as models U50N and U50CO (Figures 5, 6, 9b and 10b). Both these models exhibit radial normal regime (Aψ <60°) in the region surrounding the basin, and transition to a strike-slip regime further outward in r proj . This transition occurs at greater r proj for the model with the crustal collar (U50CO), suggesting that such loads can make important contributions to the radial nature of fault systems around SP. In the outer part of the failed zone, intermediate Aψ values lying between "pure" radial normal (30°) and strike-slip (90°) predict the presence of faults with mixed modes (Figures 9b and 10b). In this case, the mixed mode is called "trans-tension" (e.g., Dewey et al., 1998). Evidence for such a stress regime is seen in the strike-slip duplexes of the Virgil Fossae complex (Cruikshank, Umurhan et al., 2019) within the "Sputnik Radial" system. The main Virgil Fossae troughs are thought to be extensional faults (as mapped above), but the duplexes indicate right-lateral strike slip motion. At the duplex location, the strikes of the main Virgil Fossae faults veer counter-clockwise from radial, making an acute angle with the radial direction. The combination of these observations is consistent with a strike-slip stress regime with the most compressive horizontal principal stress oriented radial to the basin, exactly as predicted by the T e = 50 km models of Figures 5 and 9b, and Figures 6 and 10b.
For models with compensation of initial basin topography, we find a similarly successful subset of models: C70N ( Figure S9 in Supporting Information S1 and Figure 11b) and C70CO ( Figure S10 in Supporting Information S1 and Figure 12b) that occur at a higher value of T e (70 km) than for the uncompensated models. This finding demonstrates how basin compensation reduces the effective thickness of a shell, such that it mechanically behaves like a uniform shell of smaller thickness.
Models with shell thicknesses that differ from the favored values described above for each respective mode of compensation do not predict fault orientations consistent with those observed in the Sputnik Radial system, or have other properties that point toward their rejection as viable models. For example, basin-loading models with thinner than optimal shells (e.g., U30N, Figure S3 in Supporting Information S1, Figure 9a, and C60N, Figures 11a) demonstrate only a narrow (∼100 km width) annulus of concentric normal fault regime (Figure 9a), with a robust strike-slip regime seen for r proj ≳ 550 km. This finding is at odds with the observed radial nature of the Sputnik Radial system (Figure 3). Further, for the basin-load model U30N the vertical extent of the failed region is quite deep, constituting half of the shell thickness for 390 ≲ r proj ≲ 610 km and more than 80% of it near r proj = 430 km ( Figure S3 in Supporting Information S1). Such a finding would suggest pervasive faulting and/or deep-seated tectonic disruption in this region, inconsistent with the moderate extent and topographic expression of faulting around SP. And while adding crustal collar loads increases the amount of predicted concentric normal faulting (Figures 10a and 12a), these zones are interrupted by an annulus of ambiguous fault orientation (white in Figures 10a and 12a), and the depth of the predicted failed zone in the uncompensated model U30CO is even greater in depth and width ( Figure S4 in Supporting Information S1) than the model of Figure S3 in Supporting Information S1, compounding the misfits.
Models with thicker than optimal shells also fail to match observed fault distributions in the "Sputnik Radial" system. Consider the basin-load models U70N ( Figure S5 in Supporting Information S1 and Figure 9c) and C80N (Figures 11c): there are two problems. First, substantial fractions of the predicted failed zones are occupied by a concentric normal faulting regime (Aψ <0°, orange-red to red colors), whereas no faults with such orientations are seen in the Sputnik Radial system (or elsewhere around SP). Second, the relatively narrow (∼150 km, Figure 9c) and extremely shallow (around 1 km, Figure S5 in Supporting Information S1) extents of the predicted failure zones for the models lacking a crustal collar load are inconsistent with the observed radial ranges and robustness and topographic expression of faulting. The addition of a crustal collar ameliorates these problems, but only to a small degree (models U70CO, Figure S6 in Supporting Information S1, Figure 10c, and C80CO, Figure 12c). Similar arguments apply to the "West Sputnik" system (red in Figure 2). The orientation relative to SP varies considerably along its length: where the system intersects the edge of SP near 37°N 145°E in northern Viking Terra, it shows predominantly radial orientations from about r proj = 550-700 km, but as it progresses to the southwest, the fault strikes diverge to an oblique angle of 25° or so counter-clockwise from radial in Cthulhu Macula. As stated above, this divergence is consistent with mixed-mode trans-tensional faulting (as indicated by the Aψ values in Figures 5b, 6b, 9b and 10b), but also may reflect contributions from global-scale stress sources such as TPW (e.g., Keane et al., 2016). Overall, uncompensated models with T e = 50 km and compensated models with T e = 70 km provide the best fits to the orientations of the West Sputnik system.
The more diffuse Northwest Sputnik and East Sputnik systems (yellow and dark blue lines in Figure 2, respectively) generally have shorter fault segments than the West Sputnik and Sputnik-Radial systems. Their fault strikes are generally quite consistent across their extents, which means that both feature zones with fault strikes that are to radial to SP (the northernmost faults in Northwest Sputnik and the southernmost ones in East Sputnik) and zones where strikes are oblique to SP (the southernmost faults Northwest Sputnik and the northernmost ones in East Sputnik). The components of these fault systems that exhibit radial orientations are consistent with the predictions of the uncompensated models U50N (Figures 5 and 9b) and U50CL (Figures 6 and 10b) and the compensated models C70N ( Figure S5 in Supporting Information S1 and Figure 11b) and C70CO ( Figure S6 in Supporting Information S1 and Figure 12b). Components that exhibit oblique orientations indicate different conditions than for the optimal-thickness models, which could include: (1) a stress field that favors formation of obliquely oriented strike-slip faults; (2) a stress field that favors formation of concentrically oriented normal faults; and (3) inherited fabric from pre-basin stress state or fault system, and either faults were reactivated by the Sputnik basin loading stresses or new faults formed by the superposed effects of the old and new stress fields. Option 1 is consistent with the results of models such as U30N ( Figure S3 in Supporting Information S1 and Figure 9a) and C60N (Figures 11a) in which the circum-basin failure zones are dominated by a strike-slip stress regime, suggesting that lateral variations in lithospheric thickness around the Sputnik Basin could play a role in the observed tectonic signatures. Option 2 would be consistent with the results of models U70N ( Figure S5 in Supporting Information S1 and Figure 9c) and C80N (Figures 11c) in which stress regimes favoring concentric normal faults occur at the basin margin, again suggesting heterogeneity in shell thickness. However, the continuity of trends between the proximal and distal faults in the northern part of the East Sputnik system (Figures 2  and 3), for example, suggests a common stress field for both and not a transition from concentric normal to radial normal to strike slip as predicted in the models of Figures 9c and 11c, thus casting doubt on this option. Option 3 is appealing given the broad consistent orientation trends within the Northwest Sputnik and East Sputnik systems (faults in both systems are primarily oriented NW-SE), and the apparent continuity of the latter with faults in the basin-proximal West Sputnik system (Figures 2 and 3).
We do not consider the relationship of the RTS to SP loading because of evidence that this system predates the SP impact . However, it is clear that the Sputnik basin has interacted with the older RTS where the disrupted blocks of material (such as al-Idrisi Montes) are seen to break off the water ice shell and embed within the nitrogen ice (e.g., Moore et al., 2016;O'Hara and Dombard, 2021;Skjetne et al., 2021;White et al., 2017).

Interplay of Flexural and Membrane Support
The stress and deformation of Pluto's icy shell from the nitrogen ice loading of the Sputnik basin, which are the important properties of the modeled stress states presented above, result from interplay between two important modes of lithospheric loading: a bending, or "flexure" response, and a stretching, or "membrane" response (Banerdt et al., 1992;Janes and Melosh, 1990;Turcotte et al., 1981). In general, a membrane-type response becomes more prominent with increasing contribution of sphericity to the response (e.g., Turcotte et al., 1981), a function of the ratio of characteristic load width to planetary radius (which can also be expressed as degrees of arc). However, some amount of stretching response can be seen in planar-geometry flexure solutions as well, particularly for more physically complete "thick plate" formulations (e.g., Comer, 1983). Nonetheless, membrane stress components induced by Sputnik basin-sized loads can exert a primary control on the tectonic expression of a planet, primarily through the extensional contribution to the out-of-plane stress (see below).
Both flexure and membrane responses are considered "tectonic" stress components as defined above, because they are superposed on the general state of stress induced by gravity, which we assume to be lithostatic here Bending of an elastic plate produces a characteristic "dipole" tectonic stress response, with maximum values of horizontal tectonic stress at the top and bottom (of opposite sign), with a "neutral plane" of zero stress in the middle. The dipole is the manifestation of a vertical gradient in horizontal normal stress in the absence of any stretching response. In contrast, the tectonic stresses induced by a stretching response are constant through the thickness of the lithosphere. This interplay is demonstrated in Figure 13b: the deviations from vertical (or, gradients) in the tectonic stress components represent flexural contributions, and the positions with respect to the zero-stress coordinate reflect the membrane contributions. Note that the horizontal in-plane stress components (σ Th , red lines) tend to have the larger gradients, and are thus more influenced by the bending response, whereas the out-of-plane components (σ Tϕ , blue lines) tend to be vertically more flat, and are therefore more influenced by membrane support. Further, note that the entireties of the σ Tϕ components in. Figure 13 are on the extensional (positive) side. This phenomenon is a main driver of both the radial normal fault regime zones seen in Figures 5b, 6b, 9b and 10b, and the facilitation of radial orientations for dikes beyond the basin, allowing the enhanced magma ascent potential seen in Figures 5c and 6c to be realized.
Many planetary loads that cause downward flexure are of modest lateral extent relative to the planetary radius: for example, volcanic edifices on large planets such as Earth and Venus (e.g., McGovern, 2007;McGovern and Solomon, 1998) and volcanic fill of impact basins on the Moon (e.g., Mare Serenitatis, in Solomon and Head, 1979). For such loads, the in-plane horizontal normal stress (called σ h here) is the most tensile one, leading to a favored concentric orientation of extensional faults (graben). In contrast, stretching of a spherical shell produces components of horizontal normal stress that are constant with vertical position in the shell. For sufficiently wide downward loads on elastic spherical shells, the out-of-plane normal stress σ ϕ is the most extensional, leading to a predicted radial orientation of graben (e.g., Banerdt et al., 1992;Sleep and Phillips, 1985;Turcotte et al., 1981).
Variations in model parameters can affect the balance of the flexural (bending) and membrane (stretching) responses: in general, the response is governed by whichever mechanism provides the most resistance to deformation. We find that low values of T e and high values of characteristic load width favor the membrane response, and high values of T e and low values of characteristic load width favor the flexure response. The analysis of Turcotte  (1981) reveals the physical basis for these results: the parameter σ that measures the resistance to bending is proportional to T e to the third power, whereas the parameter τ that measures the shell's resistance to deformation if bending is neglected (i.e., the membrane response) is merely proportional to T e . Thus, the ratio σ/τ is proportional to T e 2 , demonstrating that increasing T e will greatly increase the resistance to bending relative to the resistance to membrane deformation, leading to an increasingly bending-dominated response (being the stronger mechanism). Comparing Figure S3 in Supporting Information S1 (model U30N) and Figure S5 in Supporting Information S1 (model U70N) demonstrates this effect: the distal shell of the former is nearly uniformly characterized by the strike-slip regime generated by the extensional σ ϕ characteristic of the membrane response. We also note an apparent paradox: while low lithospheric thickness is stated above to favor a membrane response, the low-thickness T e = model U30N case ( Figure S3 in Supporting Information S1) exhibits high stress gradients (Figure 13b), a characteristic linked to flexural response. Conversely, the high-thickness T e = 70 km case ( Figure S5 in Supporting Information S1), favoring a flexural response, exhibits low stress gradients, a characteristic of a membrane response. These findings are in fact a natural consequence of the amount of support provided by each response under situations where both contribute: flexural deformation and stress gradients increase with decreasing T e , so the T e = 30 km case will have the highest gradients (Figure 13b), although the large magnitude of the offset toward extension of the σ Τϕ curve shows the dominant influence of the membrane response in setting the tectonic state. Conversely, the T e = 70 km case shows the lowest gradients (Figure 13b), reflecting a high level of support, but the greatly reduced magnitudes of the lateral offsets relative to the other thickness cases show the diminishing influence of the membrane response.
The combination of lithospheric bending and stretching phenomena creates the distinctive features of our models. For example, the crescent-shaped upper lithosphere fault regime zones seen in Figures 5 and 6 and Figures S3-S10 in Supporting Information S1 result from the combination of a flexural stress gradient (resulting in differing stress conditions at the top and bottom of the lithosphere) with extensional components of membrane stress that vary in such a way to (generally) yield normal faulting regimes near the margin of the basin and a strike-slip regime in the surrounding regions. However, interactions between the lithospheric thickness and the shape of the load(s) also play roles here. For example, the almost complete lack of normal faulting regime zones in the T e = 30 km models ( Figures S3 and S4 in Supporting Information S1) shows the dominant membrane nature of the response for thin lithospheres; only for a limited range of r proj does σ h become extensional in the uppermost several km of the shell (joining the higher σ ϕ ), a requirement to depart the strike-slip regime.
The work of Janes and Melosh (1990) illuminates considerations of the relative contributions of flexural and membrane loading to planetary tectonics. Consider the first sentence of their abstract: "The tectonic response of a planet's lithosphere to an imposed load contains information on the mechanical properties of the lithosphere, particularly its thickness …," a statement that encapsulates the approach we have taken in this manuscript. Janes and Melosh (1990) calculated tectonic states from surface stress components using a general model of spherical shell loading (i.e., containing both flexural and membrane components) capable of handling arbitrary shell thicknesses. These workers characterized the surface tectonic signatures of end member flexure and membrane responses, and also a "transition" response between them, for typical large (Mars) and small (Uranus's moon Miranda) planets that exhibit differing values of a support parameter q. The far-field response (radial normal fault regime) to a Sputnik Planitia-filling load at T e = 50 km ( Figure 5) corresponds to the prediction of the Janes and Melosh (1990) "transition" response for a small planet with intermediate thickness lithosphere (see the upper right portion of Figure 8b of that paper). For a reduced lithospheric thickness (T e = 30 km, Figure S3 in Supporting Information S1), the trans-basin response is instead dominated by a strike-slip regime, corresponding to the "transition" prediction for a thinner lithosphere in Figure 8b of Janes and Melosh (1990). Note that the appearance of a small and shallow radial normal regime zone in our T e = 30 km model is a result of a small positive excursion of σ r between r proj ≈ 400 and 560 km ( Figure S3a in Supporting Information S1); such an excursion is a departure from the always-negative σ r "membrane" loading state, and likely results from a combination of differences in load shape and support parameter q between the respective models.

Elastic Thickness and Compensation State of Pluto's Ice Shell With Implications for Planetary Thermal History
For models with uncompensated basin topography, elastic shell thicknesses T e significantly different from 50 km produce predictions of tectonic state that differ significantly from the observed predominantly radial normal faulting state (e.g., Figures S3-S6 in Supporting Information S1; Figures 9 and 10). The best fit estimate for such shell configurations is therefore T e = 50 ± 10 km. However, introducing compensation for the basin topography thins the shell beneath the load, thereby reducing the ability of the shell to support the load. As a result, the response at a given T e resembles that of uncompensated models at a lower value of T e (e.g., Figures 8, 11 and 12; Figures S7-S10 in Supporting Information S1). For example, the surface stress state of the compensated model C70N (Figure 9b)  in Section 6.1. Thus, the preferred range for compensated configurations is T e = 70 ± 5 km. While in a strict sense these T e ranges for uncompensated and compensated situations don't overlap, we consider these families of models to be end-members bounding the range of potential SP basin compensation states. Thus, we conclude that the thickness of the elastic part of Pluto's ice shell T e during the time that the bulk of the nitrogen ice loading occurred was between 40 and 75 km, encompassing the estimate of 60 km from limb profile topography of Conrad et al. (2021). This correspondence suggests that the time periods over which the topographic features observed on the limbs formed and the nitrogen ice load was emplaced overlapped to some extent, but does not require overlap if a quasi-steady-state shell thickness was reached at some point in Pluto's evolution. Further, this range is consistent with the evaluation of T e via heat flux considerations under expected levels of radiogenic heating of Pluto and the assumption of low shell conductivity corresponding to porous ice (Figure 3 of Nimmo & McKinnon, 2021).
Global climate modeling of Pluto (Bertrand et al., 2018) indicates that the massive scale of the Sputnik basin would cause it to be filled by nitrogen ice migrating into it from across Pluto after its excavation over a timescale of some tens of millions of years (less than 1% of Pluto's total age). If the basin is filled so rapidly, it follows that the formation of the Sputnik basin and its subsequent infilling (i.e., loading) are essentially coincidental to within some tens of millions of years. Since our T e findings correspond to the time of loading, we can then use the modeled evolution of the thickness of the elastic layer by Bierson et al. (2020) to precisely constrain the timing of the Sputnik impact using our results for the elastic layer thickness during nitrogen ice loading. According to the "hot start" model in Figure 1 of Bierson et al. (2020) (which they determined to be a more plausible model of Pluto's thermal evolution than the "cold start" model), the thickness of the elastic layer expanded from 40 to 75 km between 4.45 and 4.33 Ga, which causes us to conclude that the Sputnik basin formed at some point within this period. However, this finding assumes that the entire nitrogen ice inventory of the Sputnik basin today was available on the surface for infilling immediately following the Sputnik impact, and not (for instance) released gradually from Pluto's subsurface over a protracted timescale, which would have correspondingly lengthened the timescale of infilling of the basin.

Sub-Ice-Load Stress State
Model stress states predict thrust faulting of ambiguous orientation in the interior of the load region (Figures 5, 6 and 8; Figures S3-S10 in Supporting Information S1). However, evidence for compressional tectonics per se is lacking at the surface of the materials filling Sputnik basin. The morphology of the nitrogen ice materials is controlled by solid-state flow at current Pluto surface conditions Trowbridge et al., 2016;Umurhan et al., 2017;Wei et al., 2018) that will tend to remove tectonic evidence of lithospheric stresses, in contrast with, say, basin-filling lunar mare basalts, which do not flow at lunar surface conditions and therefore display tectonic features consistent with predicted lithospheric stress states (e.g., Solomon & Head, 1980). This is similar to the effect that detached basal boundaries have on volcano tectonics (e.g., McGovern & Morgan, 2009;McGovern & Solomon, 1993, in that transmission of horizontal compressive flexural stresses from the lithosphere into the load is inhibited by a slip boundary condition between lithosphere and load, but with the added difference that the nitrogen ice load itself removes internal differential stresses via ductile flow.

Topographic Constraints: Rim
The Sputnik basin exhibits prominent rim topography that is particularly well defined in the west profiles (Figures 4b and 14). As noted in Section 1, impact basin rim topography can arise from the response to the impact process, including structural uplift and overturn of the target rock and emplacement of ejecta at the surface (e.g., Melosh, 1989). The extent to which post-impact lithospheric loading plays a role in generating a rim signature can be evaluated by the results of our models. Models in which the sole loading is within the basin (e.g., Figure 5; Figures S3, S5, S7, S8 and S9 in Supporting Information S1) do not generate a raised rim (blue dashed line in Figure 14). Note that this result is counter to the prediction of substantial topographic "arches" generated by flexure models that neglect the effects of lithospheric curvature (e.g., Mills & Montési, 2019). The difference results from the relatively low amount of support from flexural bending, the main arch-producing mechanism (and the origin of the term "flexural arch"), for loads of SP's dimensions on a planet of Pluto's size. This effect has been demonstrated in loading models for mare-filled impact basins on the Moon, where membrane-dominated lithospheric responses produce comparably small amounts of topographic relief beyond the margin of the load: for example, Figure 9 of Solomon andHead (1979) andFigures 11-18 of Solomon andHead (1980). However, uplift resulting from a crustal collar load can produce topography at the basin margin that at least partially resembles the observed topographic profiles (red, yellow, and purple lines in Figure 14). For elastic thicknesses at the middle of our preferred range for uncompensated models (T e = 50 km), uplift from the crustal collar load produces a broad swell spanning from the basin margin to r proj ≈ 900 km. These model swells overlap with two prominent topographic highs in the west basin profiles (Figures 4b and 14): the aforementioned rim (with the modeled rim crests falling within a few tens of km of the observed rim crest at r proj ≈ 625 km) and a second outer high located between about r proj ≈ 750 and 830 km. The model profiles could plausibly contribute to the relief seen at these highs, and the outward-facing slope of the outer high is echoed in the model profiles. Of course, the prominent low between these highs disrupts the fits. However, this low likely reflects another critical aspect of the impact process: an outer basin ring (e.g., McKinnon et al., 2017;Melosh, 1989;Spudis, 1993), marked by the highest point of the outer high. Thus, it is reasonable to evaluate the fit at only the topographic highs. We conclude that uplift from a crustal collar load can contribute to the observed topographic profiles of the Sputnik basin and its surroundings, although this is not a required element given the other potential relief-generating mechanisms, including ejecta emplacement  and outward collapse of the transient crater (Potter et al., 2013).
Since these processes are not simulated in our modeling, we infer that the modeled profiles that closely match the observed rim topography (i.e., those for the nominal crustal collar loads in Figure 14) suggest that the present day rim topography is essentially a signature of uplift from a crustal collar load, with minimal contribution from ejecta to the observed topography. For the narrow and absent crustal collar load models, the profiles of which less closely match the observed topography, ejecta likely forms a more significant component of rim relief. The degree to which crater material and ejecta presently contribute to the structure of the basin and the topography of its rim is dependent not only on how the circumstances of the impact dictated the configuration of the basin structure and the distribution of its ejecta (e.g., Potter et al., 2013), but also on post-impact modification: given the evidence for widespread glacial erosion having affected the SP basin rim White et al., 2019), this process has likely contributed to significant removal of material from the basin rim and reduction in its relief.
The response of the lithosphere exerts strong controls on the possible range of crustal collar topographic expressions. The characteristic width of the uplift profile is not substantially affected by the width of the actual crustal collar load: compare the results for the nominal width (brown dashed line) and narrow (half the nominal width, yellow dashed line) crustal collar loads at T e = 50 km in Figure 14, for which the main difference is the lower height reached by the narrower (and therefore weaker) load. Thus, we conclude that the uplift wavelength is set by the response of the lithosphere, which effectively filters the wavelength of the load: thinner lithospheres will produce shorter uplift wavelengths. Also note that the best overall fit to the main rim and outer high (including outward-facing slope) comes from model U30CO (T e = 30 km), with its shorter wavelength response. However, this best fit would presume that no other process contributed to positive rim topography at the Sputnik basin, such as the aforementioned ejecta emplacement, and in any event this T e value can be rejected for the other reasons outlined above.

Topographic Constraints: Initial Basin Depth
The initial depth of the basin is also constrained by the topography. The central depth of fill d final required to match the observed −2 km elevation offset increases non-linearly with increasing initial basin depth d bc , such that merely doubling d bc for an uncompensated T e = 50 km model quadruples d final . This strong non-linearity is a consequence of the compensating lithospheric deformation resulting from a given amount of applied load: adding a load thickness increment increases the resulting surface elevation of the load by only a fraction of this increment, determined by the interplay of flexure and membrane responses as described above. On Pluto, this phenomenon is exacerbated by the extremely low compensating density contrast between shell and ocean that results in high compensating deflections relative to, say, silicate-dominated planets with large crust-mantle density contrasts. Such "diminished returns" of resultant topography from applied thickness increments means that a maximum nitrogen ice thickness reaching several tens of km is required for basins with d bc greater than 4 km. Such values of nitrogen ice thickness are unlikely on several grounds, not least that this depth is much deeper than any likely impact basin, especially as the surface of SP is already at least 2-3 km below the surrounding terrain . Further, these effects are magnified for compensated basins where the central shell thickness is decreased ( Figure S11 in Supporting Information S1). Thus, we suggest that the original depth of the basin at the time of initiation of nitrogen ice infill cannot be significantly deeper than 3 km. Such a finding is consistent with the "warm" Pluto SP impact models of Johnson et al. (2016), with post-response depths of several km, and inconsistent with the "cold" Pluto models in the same work (initial depths of 10s of km) unless significant post-impact relaxation of topography occurred. If Pluto remained cold, such relaxation would be unlikely (although impact heating of the SP substrate could contribute locally). Our results are also therefore consistent with the "hot" start and early ocean formation for Pluto found by Bierson et al. (2020).
Our models also show another primary reason why the initial shape of the Sputnik basin must be pan-shaped and shallow (less than 4 km): the amount of fill required to meet the surface topography constraint is so large that the lithosphere would be pervasively faulted through its depth in a dominantly strike-slip mode ( Figure S7 in Supporting Information S1), contradicting the observed limited amount of faulting and its dominant radial extension nature. The problem is exacerbated by the added thinning of the central shell that results from deeper basins and compensation of basin topography, yielding more deformation and stress for a given load amount. Further, possible stress/deformation ameliorating factors such as thicker lithosphere (Figures S5 and S6 in Supporting Information S1) are not viable because models with those properties fail to match observed faulting geometries.
The elevation offset of −2 km also provides a constraint on the minimum initial depth of the basin. For a panshaped basin with d bc <2.2 km, there is no solution that allows a −2 km offset, as the final relief built up within the basin will always exceed that topographic level. In the limit of no initial basin topography (the favored hypothesis of Hamilton et al., 2016), only positive relief features (essentially mountains) can result from emplacement of material with the density of nitrogen ice. In order for this scenario to match the observations, several (perhaps ranging into the high single digits) kilometers of nitrogen ice would have to be removed while more or less maintaining the underlying lithospheric deflection profile caused by the original load. Since the fluid providing the basal restoring force (the water ocean) has very low viscosity, response to load removal will be nearly instantaneous, such that a time lag between nitrogen ice load removal and basin readjustment cannot be invoked to create this situation. Thus, to create and preserve an offset, the "no initial basin" scenario requires both partial removal of the load and an increase in elastic shell thickness between times of original load emplacement and partial load removal; Hamilton et al. (2016) invoke an early epoch of "thin" lithosphere in order to facilitate creation of a deep depression from deposited nitrogen ice. Nonetheless, an "initially flat" scenario must conform to the constraints laid out above, which include a T e value at the time of nitrogen ice loading of about 50 km. The "no initial basin" scenario would therefore require a current-day T e value well in excess of 50 km, implying a high rate of internal cooling for Pluto. Further, a "no initial basin" scenario cannot explain the current elevated rim of the Sputnik basin (Figure 4), since loading will not produce such a topographic signature (Figure 7). However, the impact process itself, and subsequent deformation resulting from impact-generated loads like crustal collars, are capable of producing rim topography. Thus, the results reported here favor the interpretation of Sputnik Planitia as the site of an impact basin.
Finally, we note that the final shape of modeled basins tends toward a bowl shape, regardless of whether the initial shape was bowl-like or pan-like. For initially pan-shaped basins, the downward deflection of the lithosphere tends to remove the flat central topographic signature; compare the initial and final configurations of the nominal pan-shaped basin model (red lines) in Fig. 7c. This is a quite simple consequence of the long-wavelength natures of both flexure and membrane responses, with the greatest deflections near the center of the load and decreasing outward in r proj . Any initially flat surface will inevitably be tilted inward toward the center of the basin by such a response, and a curved shape characteristic of loading-induced subsidence will be reflected in the final basin surface configuration.

Cryomagmatism
In general, stress gradients for models with relatively thick lithospheres are too low to provide any significant enhancement of cryomagma ascent in the regions surrounding the basin. Model U70N ( Figure S5 in Supporting Information S1) reflects such a condition. However, thick-shell models with mitigating conditions such as a crustal collar load (e.g., U70CO, Figure S6 in Supporting Information S1) or basin topography compensation (C70N, Figure S9 in Supporting Information S1), or both (C70CO, Figure S10 in Supporting Information S1) offer modest enhancements, with values of Δρ efb in the 40-70 kg/m 3 range. Note that compensation can result in a inward shift of the position of the zones with enhanced Δρ efb (e.g., model C70N, Figure S9 in Supporting Information S1) to locations where cryovolcanism is not observed at the equivalent radial positions at SP (e.g., at the basin margin, Figure S9 in Supporting Information S1), although adding crustal collar loads can restore these zones to their expected locations (e.g., model C70CO, Figure S10 in Supporting Information S1).
Thus, if lithospheric stresses have aided cryomagma ascent around SP, we can infer that Pluto's T e is less than about 60 km for uncompensated shells, or 70 km for compensated ones, values also consistent with the dominantly radial extensional tectonic modes surrounding SP (see the "elastic thickness" section above).
As described in Section 1, proposed sites of cryomagmatism include Hekla Cavus, Viking Terra, Pioneer Terra, Virgil Fossae, Wright Mons, and Piccard Mons (at distances 527, 582, 815, 900, 957, and 1,254 km from the center of SP at 25°N, 175°E, respectively). The first three fall into radial ranges that overlap with zones of enhanced magma ascent in one or more or our models (Figures 5 and 6; Figures S3, S4, S6, S8 and S10 in Supporting Information S1), but the two Montes and Virgil Fossae fall beyond the zones with peak enhancement (Δρ efb > 80 kg/m 3 ). In the cases of Wright and Piccard Montes, they are closer to the south-southeast extension of SP, a zone of loading that might become a significant or even primary influence in that region, suggesting the need for further non-axisymmetric models. Virgil Fossae form part of a fault system oriented radially to SP, suggesting that the observed cyromagmatism could reflect down-strike transport in an underlying dike connected to the enhanced zone.

Basaltic Planet Analogs for Basin-Flanking Volcanism
Here, we have proposed that the emplacement of cryovolcanic centers in the region surrounding the Sputnik basin has been facilitated by lithospheric stresses induced by the basin-filling nitrogen ice load. However, associations of volcanic provinces with large impact basins have long been noted on the silicate planets Mars and the Earth's Moon, and McGovern (2018) proposed that similar stress-based magma ascent enhancement has occurred at these locations. Specific examples of such provinces include the Circum-Hellas Volcanic Province (CHVP) on Mars (Williams et al., 2009) and Mare Tranquillitatis, located within the triangle formed by the Crisium, Serenitatis, and Nectaris basins on the Moon (Litherland & McGovern, 2009;McGovern & Litherland, 2011;McGovern et al., 2014). Preservation of such provinces to the present day requires that the surface signatures of large basins and their surroundings not be removed by vigorous resurfacing that evidently characterized the larger silicate planets Earth and Venus, thereby limiting such evidence to the smaller planets, sometimes termed "one-plate planets" (Solomon, 1978). The Kuiper Belt object Pluto evidently shares such preservation of the early solar system impact record.

Future Directions
Despite the progress on understanding Sputnik Planitia and its surroundings presented here, there remain ample opportunities for improvement. Future work on the Sputnik basin can address some of the limitations of the models presented above. For example, the significant departure of the SP nitrogen ice unit from axial symmetry in the south-southeastern extension of the main basin, which may reflect an elliptical basin shape or irregular low-angle impact geometry, is unaccounted for here. While we do not expect accounting for such variations to substantially change the main conclusions reached here regarding lithospheric thickness and basin dimensions, implementing models that account for non-axisymmetric shape of SP could address in particular the role of loading in enhancing cryomagmatic potential at the locations of the prominent Wright and Picard Montes. A 3-dimensional approach can account for such asymmetries, allowing for explicit calculation of tectonic and topographic expressions from the full extent of the nitrogen ice load (e.g., McGovern et al., 2020). Calculation of time-dependent responses to loading can also increase our understanding. For example, implementing viscoelastic rheology would allow the response of the whole shell to be accounted for, including the warm ductile lower regions. Further, models with plastic rheology can account for the stress relief from brittle failure (faulting), allowing for a more self-consistent evolution of the surface and interior stress states. Finally, we used a value for ocean density that did not account for salinity, which can make oceans denser (by perhaps as much as 200 kg/m 3 : see Johnson et al., 2016). A denser ocean than modeled here would result in a greater amount of buoyant support, resulting in lower deflection, stress, and stress gradient values for a given load. Such a variation might relax constraints on initial basin depth and shape and change (likely reduce) favored ranges of T e inferred from the models presented above. The extent of such differences can be tested by future models using a plausible range of ocean density values (for which our value is the lower bound).

Conclusions
The fault systems that surround the Sputnik basin are critical clues to the characteristics of the mass of nitrogen ice that fills the basin and to the properties of Pluto's icy shell lithosphere. Numerical models of nitrogen ice loading within SP predict distributions of faulting with mechanisms (dominantly normal, but with some transtensional contributions) and orientations (dominantly radial, with increasing obliquity with increasing distance from the basin center) that match the observed distributions of tectonic features surrounding SP. These distributions are strongly diagnostic of the elastic thickness of Pluto's spherical ice shell: preferred values fall in the range 40-75 km, with increasing amounts of compensating relief on the OSB favoring the higher end of that range. Models with elastic ice shell thicknesses out of the preferred range for a given compensation state predict tectonic states that differ significantly from the observed predominantly radial normal faults. The initial depth of the basin is also constrained by the topography. The required depth of fill increases non-linearly with increasing initial basin center depth d bc , such that merely doubling d bc for the T e = 50 km model sextuples the load thickness required to reach the observed −2 km elevation offset. Thus, we suggest that the original depth of the basin cannot be significantly deeper than 3 km, a finding consistent with the "warm" Pluto SP impact models of Johnson et al. (2016), and therefore also with the "hot" Pluto origin of Bierson et al. (2020). We also find that the initial basin cannot be shallower than 2.2 km; otherwise the elevation offset could not reach −2 km depth. This finding rules out scenarios for which there was no pre-existing basin (e.g., the favored scenario of Hamilton et al., 2016), barring a fortuitous combination of shell thickness evolution and load erosional history. Uplift generated by a "crustal collar" load may contribute to part of the observed topographic signal of the basin and its surroundings. The combination of extensional out-of-plane lithospheric stress and positive vertical gradients of tectonic stress generated by the loading creates an enhancement of cryomagma ascent in dikes in an annular region surrounding the basin, corresponding to the locations of many prominent sites of proposed cryovolcanism on Pluto.

Data Availability Statement
The ArcGIS shape files and auxiliary files for the eight mapped tectonic systems, the depressions of the RTS, the boundary of Sputnik Planitia, and the approximate topographic rim of the Sputnik basin (as shown in Figures 1  and 2), are included in a zipped folder included with the Supporting Information entitled "ds01.zip," available in the Figshare external data repository . The topographic profiles in Figures 4 and 14 and the COMSOL Multiphysics (Version 5.6) model results from Figures 5-13 and 14, and corresponding data and model results in the Figures S1-S11 in Supporting Information S1, have also been made available on the Figshare external data repository . The base maps we used include the LORRI-MVIC global mosaic and the global stereo DEM of Pluto, both projected at 300 m/pixel . These are archived in the PDS Imaging and Cartography Node, and can be downloaded at the following links: global mosaic: https://astrogeology.usgs. gov/search/map/Pluto/NewHorizons/Pluto_NewHorizons_Global_Mosaic_300m_Jul2017. global DEM: https:// astrogeology.usgs.gov/search/map/Pluto/NewHorizons/Pluto_NewHorizons_Global_DEM_300m_Jul2017.