Propagation of Vertical Fractures through Planetary Ice Shells: The Role of Basal Fractures at the Ice–Ocean Interface and Proximal Cracks

The presence of smooth, young surfaces indicates that regions of Enceladus and Europa have been resurfaced through recent or ongoing activity related to the eruption of liquid water from subsurface reservoirs. For interior material to erupt or flow out onto the surfaces of these satellites, fractures would have to vertically penetrate the ice shell to the depth of a subsurface reservoir or ocean. Here we use linear elastic fracture mechanics to show that accounting for fracture interactions makes it much more difficult for fractures to penetrate the entire ice shell than previous estimates. We found that fractures that originate from the surface are unlikely to penetrate the entire shell thickness, even for the upper range of tectonic stresses estimated for each moon. Tensile fractures that initiate from the bottom of the icy shell—as observed in terrestrial ice shelves—propagate further into the icy shell than surface crevasses but still do not penetrate the entire ice thickness. However, full ice shell thickness fracture is possible if shear failure connects the surface with deep-penetrating basal fractures in thinner ice shell thicknesses and under certain stress conditions. This suggests that the combination of tensile and shear failure may be important and necessary for the formation of a connection from the surface to the ocean below.


Introduction
Recent observations have revealed that liquid water may be erupting onto the surfaces of several icy satellites (e.g., Hansen et al. 2006;Porco et al. 2006;Roth et al. 2014;Sparks et al. 2016Sparks et al. , 2017Dalle Ore et al. 2019), a process that, along with putative tectonic overturn (Kattenhorn & Prockter 2014), could maintain their geologically young appearance. The clearest example of such eruptions are the jets of water escaping from the tiger stripes, a set of parallel fractures at the south pole of Saturn's moon Enceladus (Porco et al. 2006;Spencer et al. 2006). On Jupiter's moon Europa, the recent discovery of a possible water vapor plume (Roth et al. 2014;Sparks et al. 2016) emanating from the southern hemisphere, initially suggested to be associated with tensile fractures, supports the likelihood that the moon maintains active resurfacing processes. In both cases of extrusive and effusive release of water, a conduit must exist between a subsurface reservoir and the surface, regardless of whether that reservoir is a global subsurface ocean or a shallow localized water pocket. Observations show that the surfaces of the satellites are highly fractured; it may be possible that some fraction of the observed surface features reach some form of subsurface water reservoir.
Recent studies of fracture propagation in icy satellites have used linear elastic fracture mechanics (LEFM) to examine the penetration depth of surface fractures (e.g., Lee et al. 2005;Qin et al. 2007;Rudolph & Manga 2009;Walker et al. 2012;Walker & Schmidt 2015). Previous estimates of fracture propagation speeds (Rhoden et al. 2010) suggest that LEFM is appropriate under the assumption of quasi-static crack growth (Sharon & Fineberg 1999). Accounting for variable porosity, Lee et al. (2005) and Qin et al. (2007) showed that a surface fracture could penetrate a 0.33-1 km thick ice shell with tectonic stresses <200 kPa. Alternately, using a model that assumed a dual-layer ice shell to account for stress relaxation in the lower reaches of the shell, Rudolph & Manga (2009) showed that under a total stress of 3 MPa a fracture could completely penetrate a 5 km Europa ice shell and a 30 km Enceladus ice shell. Although these studies provided estimates of possible fracture depth, they focused on the propagation of an isolated fracture. Observations show dense fracturing in most icy satellite surfaces (e.g., Pappalardo et al. 2009;Spencer et al. 2009), and as such, the parameterization of single isolated fractures may be unrealistic for determining fracture propagation depth. Observations in terrestrial ice and other materials have shown that interaction between closely spaced fractures (stress shielding) must be taken into account to assess accurate penetration depths (Karihaloo 1979;Lam & Phua 1991;Zou et al. 1996; Van der Veen 1998; Thomas et al. 2017). Moreover, previous studies have focused on the penetration of surface fractures because warmer ductile ice at the base was considered too viscous to fail brittlely (e.g., Lee et al. 2005;Rudolph & Manga 2009). However, terrestrial observations (e.g., Humbert & Steinhage 2011;Luckman et al. 2012;McGrath et al. 2012a;Humbert et al. 2015;Walker et al. 2021) show intense basal crevassing and rifting in ice shelves and floating glaciers, where they initiate owing to mechanical heterogeneity (e.g., Luckman et al. 2012;Gaume et al. 2017), or longitudinal stress gradients (e.g., Anderson & Grew 1977;Maimon et al. 2012). Until the onset of convection, basal ice in planetary bodies likely has similar properties to warm ice (e.g., basal marine and meteoric ice) on Earth. Due to its small size and low gravity, the base of Enceladus's ice shell may be in an analogous state to terrestrial basal ice. For Europa, the transition to convective from conductive heat loss is expected once the ice shell grows beyond ∼10 km (e.g., McKinnon 1999). Even on Earth there is observational evidence showing that basal cracks, crevasses, and fractures form in warm base/polythermal glaciers, particularly where contaminants from the base material (e.g., muddy water, till) can disrupt the ice structure (Nye 1976;Fountain et al. 2005;Hambrey et al. 2005;Hewitt & Schoof 2017).
We conduct an LEFM analysis to determine upper bounds on through-ice shell fracture by considering the effects of fracture spacing on basal and surface cracks. We first describe our model of the propagation of surface cracks and modify the single-crack LEFM model by incorporating fracture interaction. We then discuss the formation of fractures that originate from the bottom of the shell (basal fractures) and show that it is easier for basal fractures to approach full-thickness rupture of the icy shells. Lastly, we show that the requirement for fractures to penetrate through the entire ice shell thickness may be unnecessary because shear failure along faults may enable basal fractures to connect the surface and ocean.

Linear Elastic Fracture Mechanics
Fracture mechanics is based on the assumption that all materials contain defects that affect their load-bearing capacity. LEFM is a simplified approach that models the propagation of initial starter cracks or flaws in an elastic layer by assessing stress concentration at the crack tip. Ice is not a linear elastic material; its deformation is best described by a viscoelastic rheology. However, on short timescales-like that of fracture propagation-ice behaves approximately elastically. As such, LEFM is an adequate approximation to study crack propagation in ice (e.g., Rist et al. 1999). In LEFM theory, fractures propagate when the concentration of stresses near the crack tip -the stress intensity, K I -exceeds a material parameter called the critical stress intensity or fracture toughness, K Ic . Figure 1 displays the model scenario investigated here, featuring series of both surface and basal cracks in a layer of ice. The total stress intensity at a crack tip is a function of the combined stresses acting on a crack. As such, if these stresses are known, the distribution of stress near the crack can be described. The stress intensity factor K I for Mode I (opening) fracture subjected to any applied stress is generally calculated as (e.g., Broek 1982) bs where d is the crack length ( Figure 1) and β is a dimensionless quantity that depends on crack geometry, whose solutions have been calculated for a variety of crack scenarios and published in handbooks, e.g., Tada et al. (2000). The term σ xx represents the total stress acting on the crack face to open the fracture, and the net stress intensity factor K I is obtained through superposition of contributions from each of the stresses present (Van der Veen 1998). The net stress intensity at a crack tip in an ice Figure 1. Left: cartoon illustrating the ice shell crack problem investigated here. An ice shell of thickness H is subjected to a tectonic (resistive) stress of R xx , and at depth overburden pressure P O . Surface cracks of depth d and basal cracks of height h exist in the shell, and spacing between adjacent cracks is 2W. Basal fractures are assumed to be filled with water from the ocean below, thus subjected to water pressure P W . The stress intensity factor is computed at the crack tip (dashed circles) assuming that initial voids/cracks are larger than the plastic zone, r p > 15.9 m. Any uncracked layer between the surface and basal crevasses is of thickness L S , through which shear fractures (red dashed) may connect the surface with basal crevasses of sufficient height h. Further, fracture tip stress concentration and thus propagation behavior are modified by nearby fractures; gray dashed lines show transition around the fracture tip from shielding (e.g., where same-plane fractures lie), to amplification (e.g., fractures approaching each other from opposite planes), which occurs at approximately 70°from the direction of propagation. Right: profiles of stresses at depth considered here to modulate the stress intensity at the crack tip. Depth-independent extensional tectonic stress σ T acts to open a crack, while lithostatic (overburden) stress σ L acts to close the crack increasingly with depth. If there are water-filled basal cracks, water pressure σ w exerts back pressure on crack walls, acting to oppose overburden pressure and open the crack, decreasing with height above the base.
shell is decomposed into (1) a term related to the far-field depth-independent tectonic stress in xx, σ T ; (2) the depthdependent ice overburden, σ L ; and (3) water pressure within the fracture (basal fracture case only), σ w . The latter two are often considered in models of ice shell fracture. Crack propagation occurs when the net stress intensity surpasses K Ic (i.e., when K I > K Ic ). Upon crack initiation, an estimate for the depth to which a fracture will propagate is obtained by finding the depth at which K I K IC . The critical stress intensity (K IC ) most often used for water ice is 0.15 MPa m 1/2 (Rist et al. 2002;Litwin et al. 2012). We experiment with a fracture toughness range of 0.1-0.4 MPa m 1/2 (Schulson & Duval 2009), consistent with that value and those of previous icy moon fracture studies (Lee et al. 2005;Rudolph & Manga 2009). In reality, elastic stresses do not become infinite at the crack tip because plastic deformation occurs and this maintains a bounded state of stresses. LEFM theory is applicable as long as the length of a fracture is much larger than the area in which plastic deformation occurs near the crack tip, r p . An estimate of r p , dependent on plastic yield stress σ y and stress intensity K I , can be calculated as the location at which the elastic stress equals the plastic yield stress (Broek 1982) as To illustrate that LEFM is appropriate for our problem of ice shell crack propagation, we first estimate r p in a sample ice shell in the case where K I = K IC . The yield stress of ice is a crucial piece of information but is not well constrained. Previous models, for instance, have assumed the yield stress at the surface of Europa to be just 0. The ice shells of Europa and Enceladus, among others, are generally considered to be in an extensional stress regime, at least in the near surface and in the absence of refreezing at the base (e.g., Prockter et al. 2002;Kattenhorn & Hurford 2009;Kattenhorn & Prockter 2014). Here, we consider the far-field tectonic stress (σ T ) to be tensile in nature and depth independent.
The stress intensity due to σ T in the shell, K I T ( ) , acts to open a fracture and is expressed as Here, fractional penetration depth l = d H (fracture depth d, ice thickness H). As H increases, λ → 0, and K I T ( ) would then approximate the expression by Weertman (1973) and Smith (1976Smith ( , 1978 for propagation of an edge crack in an infinite half-space (H not finite). It has been shown that the infinite half-space model underestimates the depth to which a fracture penetrates, as opposed to the nonlinear model (Equation (4)) in terrestrial glaciers (Mottram & Benn 2009).
Second, lithostatic stress, or ice overburden pressure, is a depth-dependent quantity that describes forces acting to close the fracture (Weertman 1973). It is likely that the density of ice comprising the ice shells of interest is more porous in the near surface than it is in the bottom layers. Nimmo et al. (2003) estimated that higher porosities persisted through the upper half of Europa's elastic layer before going to zero in the lower half of the ice shell. We present results using constant depthaveraged density and also for variable density to demonstrate the effects of assumed density on understanding failure modes of an ice shell.
Assuming a density profile ρ i (z), the overburden stress σ L at depth can be expressed as where z is depth below the surface s. For a constant density profile, Hartranft & Sih (1973) and Tada et al. (2000), the stress intensity factor for overburden stress, K I L ( ) , acts to close the crack, and it is expressed as . The combination of stress intensity factors due to tectonic stress (K I T ( ) ) and lithostatic stress (K I L ( ) ) are often used to determine the penetration depth of fractures. However, another factor that contributes to the stress concentration at a given crack tip is the spacing of nearby fractures. Nearby fractures, where spacing is much smaller than ice thickness, have a diminishing effect on the stress concentration at the crack tip. The stress intensity factor due to tectonic stress in multiple fractures (K I T m , ( ) ) spaced width 2W apart can be expressed as (Benthem & Koiter 1973) S is a depth-independent function of the ratio of fracture halfspace W to penetration depth d, expressed as We note that solutions in Equations (4), (8), and (10) are based on polynomial curve fitting, highly accurate for the geometry and boundary conditions used to achieve those fits (e.g., Broek 1982). However, this approach was used for thick glacier ice (100 s m) by Van der Veen (1998) and resulted in estimates that were accurate within an order of magnitude of observations. The purpose of this study is to provide bounds on through-ice shell fracture; thus, this order-of-magnitude accuracy is appropriate for our purposes. In the case of an isolated surface crack, through superposition, Equations (3), (7), and (9) allow us to determine the net stress intensity factor at the crack tip K I . We calculate the depth to which the fracture will penetrate into the ice shell by calculating the depth at which K I K IC . We note that the accuracies of these approximations are 0.5% for any λ in Equations (3) and (7) and up to 1% for any S in Equation (9) (Benthem & Koiter 1973;Tada et al. 2000).

Basal Cracks
The lower reaches of icy shells have been considered to viscously relax any stresses present aside from overburden pressure (Manga & Wang 2007;Rudolph & Manga 2009;Yin et al. 2016). Likewise, terrestrial floating ice shelves are approximated as viscoelastic solids, with warmer, more ductile ice at the base; however, basal crevasses are ubiquitous (e.g., Luckman et al. 2012;McGrath et al. 2012a;Bassis & Ma 2015;Logan et al. 2017). Assuming that basal fractures initiate upon buildup of elastic stress and propagate on an elastic timescale, we can also apply LEFM techniques to predict basal crevasse heights, ignoring effects of viscous relaxation (e.g., as in Sandwell et al. 2004;Smith-Konter & Pappalardo 2008;Craft et al. 2016).
Overlying a subsurface ocean, any bottom crack will become filled with water. Fluid within a crack will apply pressure on the walls that opposes overburden pressure (e.g., Weertman 1973). Thus, it can alter the height to which an otherwise dry fracture could penetrate. Assuming hydrostatic equilibrium, water at the base of the ice shell would rise up through any fracture to the piezometric head. The water pressure (σ w ) in a fracture at a height z above the base of the shell can be approximated by (e.g., Van der Veen 1998) where ρ w is the density of water filling the fracture, H p is the height of the piezometric head, and z = 0 at the base of the ice shell. Water pressure decreases as water rises in the shell toward H p , and it is zero above H p . At the base, the stress intensity factor in a basal crack subjected to lithostatic and water pressures can be expressed as (Van der Veen 1998) Here, G(γ,λ) is defined in Equation (8), and h is the height of the basal fracture. The importance of fracture interactions is not limited to surface fractures; the same closely spaced fracture model can be used to estimate basal fracture heights. Similar to the surface case (Section 2.2.1), through superposition the net stress intensity factor can be computed by summing the contributions from tectonic stress (K T 1 ( ) ) and lithostatic and water pressures at the base ( In the case of an isolated basal crack, we determine the height to which it will penetrate into the shell where In the multiple closely spaced fracture case at the base, we determine the height to which a single fracture could penetrate into the shell as

Shear Failure
The main goal of this study is to determine bounds on the conditions for, and potential of, tensile fracture propagation to enable ocean-surface conduits. However, we are not strictly limited to tensile failure scenarios. If tensile fractures do not propagate all the way through the shell to connect the surface and ocean, shear failure along optimally oriented faults may enable the complete failure scenario. Although shear stresses act to separate the fault planes, the normal stress acting on the fault and its frictional properties act to keep the fault closed (e.g., Byerlee 1978). Typically, the normal stress acting on a fault is the lithostatic stress, or overburden pressure, σ L , as described previously. Here, we also consider the diurnal tidal stress found at both Europa and Enceladus, where modeled values of resolved shear stress along faults suggest peak magnitudes between τ s = 80 and 100 kPa, respectively (Hoppa et al. 1999;Nimmo et al. 2007;Smith-Konter & Pappalardo 2008;Rhoden et al. 2010). Given the increase with depth of the overburden pressure, shear stresses acting to separate the fault plane would likely only exceed the normal stress at relatively shallow depths. As this study is mainly focused on the role of basal fractures in ice shell dynamics, we consider an intact layer of ice of thickness L S above the water-filled, stress-free crack to be the most likely region for such shear slip to allow full ice shell failure through a combination of elastic and shear failure (e.g., Figure 1). Considering tidal stresses that are active on both Enceladus and Europa, we can determine the stress magnitude at which complete failure of layer L S would occur by determining when the depth-integrated tidal stress s t exceeds the depth-integrated differential stress t D . This thickness of ice L S must support the same depthintegrated tidal stress s t as the global tidal stress σ t acting over the uncracked thickness H. That is, and as such the stress concentrated within the layer L S is an amplified version of the regional tidal stress acting over intact segments of the ice shell elsewhere (Smith-Konter & Pappalardo 2008;). As in Turcotte & Schubert (2002), the stress acting to resist crack opening can be expressed as a function of overburden and the coefficient of friction Integrating over the layer of thickness L S , Using Equations (14) and (16), we can determine a bound on the stability of an ice shell against complete failure by determining the upper limit for the depth to the top of a basal fracture L S that could be stable against shear failure as the point where the depth-integrated tidal stress s t exceeds the depthintegrated yield strength t D , or up to the point where the remaining intact layer over the basal fracture follows as which is similar to the expression for the upper bound on crack depth determined by Sandwell et al. (2004) in terms of tidal strain. This suggests that the portion of the ice shell L S = H-h must be less than the value on the right-hand side of Equation (17) to allow for complete ice shell failure. It is possible that an optimally located surface fracture exists in the vicinity of the fault, which would effectively decrease the depth L S ; in this case, Equation (17) still serves an upper bound on the required thickness for stability. This approach ignores the potential for interactions between surface and bottom fractures, as most interactions will diminish over time as ductile flow of the ice relaxes stress near crack tips. Hence, our expression for an upper bound L S is appropriate so long as the timescale for fracture emplacement is long compared to the timescale over which ductile flow is important. Over shorter timescales, we must consider interaction effects (Thomas et al. 2017). In that case, Tada et al. (2000) and Peng & Jones (2015), among others, show that coplanar fractures propagating toward each other are either amplified or shielded from the stress concentration around the other. Gray dotted lines in Figure 1 show zones of amplification and shielding ahead of a crack tip; this transition has been experimentally derived to be at approximately θ = 70° (Gong & Horii 1989). Following the methods of Horii & Nemat-Nasser (1985) and Gong & Horii (1989), we can estimate the amplification of K I as the fractures approach each other as q q q q q q = + + - 11 cos 8 cos 2 3 cos 3 4 sin 2 6 sin 3 2 10 sin 5 2 . 18 This expansion is based on curve fitting; this zeroth-order case assumes that cracks are similar in scale. Parameter c is the half-length of the opposing crack, and d is the offset distance. Despite this magnifying effect, we find in our cases presented below that through-shell tensile fractures remain unlikely at reasonable stress values. As such, we are able to assume that our model is appropriate given either that the interaction of surface-basal fractures occurs over a long-enough timescale that such stress relaxes or that the amplification factor over short timescales (Equation (18)) does not significantly affect depths.

Model Parameters and Experimental Setup
Determining penetration depth primarily requires an initial assumption of ice shell thickness, density, and the magnitude of the tectonic stress. There are a variety of sources of tectonic stress at icy moons that may contribute to the propagation of fractures. These include the diurnally varying tidal stresses, which, due to a Laplace-resonant orbit and orbital eccentricity at Europa and Enceladus, respectively, are on the order of 0.1 MPa (e.g., Hurford et al. 2007Hurford et al. , 2009). The cooling of icy material at the base of the ice shell may induce extensional stresses in the upper part of the shell that are of order 10 MPa in shells >10 km thick (Nimmo 2004). These stresses, however, are likely reduced to less than 3 MPa for Europa and less than 10 MPa for Enceladus when compressibility of the subsurface ocean is considered (Manga & Wang 2007). Europa may have experienced nonsynchronous rotation (e.g., Leith & McKinnon 1996) or undergone true polar wander (Schenk et al. 2008). In these cases, induced tensional stresses would be of order 1-10 MPa. It has been suggested that diapir-induced reorientation of Enceladus's ice shell may have induced tensile stresses of order 10 MPa. Additionally, Patthoff & Kattenhorn (2011) highlighted features at the south pole that suggest nonsynchronous rotation of the shell, with induced stresses between 1 and 5 MPa. Our LEFM approach, by definition, is a static fracture analysis method, meaning that the time-evolving nature of these applied stresses is not considered. In particular, as discussed in Rudolph & Manga (2009), among others, there is a large portion of the ice shell that is capable of stress relaxation, which would not allow for stresses to accumulate instantaneously. This fact has the effect of lowering the magnitude of the maximum stress expected. Here, we use the total stress values noted above with the intent of demonstrating upper bounds on ice shell fracture penetration, but we stress that these are not the actual stress magnitudes likely to accumulate instantaneously in the ice shell. Thus, the results of our analysis presented here signify upper bounds on fracture penetration through different ice shell thicknesses under a range of stresses that either currently act or at some point may have acted on the ice shell.
Estimates for the thickness of Europa's ice shell range between 1 and 30 km (e.g., Nimmo et al. 2003;Billings & Kattenhorn 2005;Lee et al. 2005;Rudolph & Manga 2009). Estimates for Enceladus's ice shell thickness vary, when considering only the south pole or the ice shell on the whole, between 5 and 90 km (Schubert et al. 2007;Patthoff & Kattenhorn 2011;Beuthe et al. 2016;Čadek et al. 2016;Le Gall et al. 2017). Thinner estimates are favored particularly in models of the formation of the south polar terrain and the tiger stripe fractures (Nimmo et al. 2007;Smith-Konter & Pappalardo 2008;Olgin et al. 2011;Walker et al. 2012).
In our models, we considered H = 5-30 km for Europa and H = 5-90 km for Enceladus to represent bounds on shell thickness estimates in the literature. We test a range of stress values likely previously or currently active on the satellites, encompassing the range of published estimates. We assume that porosity changes as a function of depth as in Johnson et al. (2017) and references therein. As an upper bound, we considered surface layer porosities of 33% for Europa (Black et al. 2001;Lee et al. 2005) and 20% for Enceladus (e.g., Kieffer et al. 2006), both of which go to zero in the lower parts of the ice shell (e.g., Johnson et al. 2017).
First, we consider an isolated fracture (surface or basal) in each ice shell and determine the fractional penetration depth (d/H) under this range of conditions. Next, we model the same fracture under the same range of conditions, but instead set it among adjacent fractures at a range of separation distances. We varied these separation distances between 10% and 50% of the total ice shell thickness, based on observations of ice shelf basal crevasses in terrestrial ice shelves  Leonard et al. 2017;Hemingway et al. 2020). We then follow these results with an analysis of how shear failure might enable surfaceocean conduits in the case that they are not predicted by tensile failure alone.

Results
Using the model setup described above, we determine variations in fracture behavior in the ice shells of Enceladus and Europa. First, we determine the conditions under which surface and/or basal fractures may completely penetrate a given thickness of ice. Here, we define "complete fracture" as d/H (or h/H) >95%, informed by error margins. Second, in the case that tensile fractures fail to do so, we determine the conditions under which shear failure along a fault in the remaining unfractured ice could create a connection. Zoomed-in insets demonstrate that (a) the difference in depth estimates could exceed ∼250 m owing to a change in density profile in a 5 km Europa ice shell, while (b) in a 5 km Enceladus ice shell estimates could be offset up to ∼150 m. Reduced shallow density allows this slightly further penetration into the shell owing to the decreased overburden pressure; the effect is more pronounced near the middle of the ice shell (sections highlighted in (b) and (c)).

Isolated Surface Fractures
To calculate the depth to which surface fractures penetrate, we determine the location at which Figure 2 not only highlights the effect of assuming constant versus depth-dependent density, showing that for a lowerdensity upper layer there is a slight increase in propagation depth, but also depicts the effects of acceleration due to gravity. Enceladus's lower gravity, g = 0.113 m s −2 , is an order of magnitude less than Europa's, g = 1.315 m s −2 , and allows for increased propagation depth since the lithostatic overburden is less. Figure 3 shows fractional penetration depths of a single surface crack on Enceladus and Europa under a range of tectonic stresses and a reduced upper shell density. Errors associated with the polynomial curve fitting for this application (Section 2.2.1) are reflected by the shaded regions, which illustrate the resultant depths across a range of applied stress and ice shell thicknesses. An Enceladus ice shell would have to be less than 30 km thick and subjected to at least ∼2-4 MPa of tectonic stress to completely fail (d/H > 95%) in the case of a single fracture. At Europa, an ice shell less than 5 km thick could completely fracture (d/H > 95%) under ∼2-3 MPa, consistent with results of Rudolph & Manga (2009). Even in the thinnest ice shells considered, tidal stresses (∼100 kPa) are not enough to propagate surface fractures through the entire ice shell thickness. A single surface fracture in a 5 km thick Enceladus ice shell would propagate roughly ∼2 km deep. Similarly, in a 5 km Europa ice shell, an isolated surface fracture would propagate 100 m deep.

Multiple Adjacent Surface Fractures
In the case of adjacent fractures (nonisolated), the magnitude of stress required to propagate such fractures is significantly larger. This is demonstrated in Figure 4, which illustrates the fractional penetration depths for a variety of fracture spacings between 10% and 50% of the ice shell thickness. In each case, decreased fractional penetration occurs when fractures are closer together. Comparing results overall (Figure 3 vs. Figure 4), significantly less penetration occurs in fractures set among an array of fractures than in the single fracture case. For example, the minimum fracture spacing tested was 0.10H. In an H = 5 km Europa ice shell (spacing 2W = 0.5 km), under 10 MPa stress a fracture would penetrate roughly 43% of the thickness, or ∼2.1 km deep. The maximum fracture separation tested was 0.5H; in the 5 km ice shell (spacing 2W = 2.5 km), under 10 MPa stress a fracture in this system would penetrate ∼60% of the thickness, or ∼3 km. We note here that much of Europa has been imaged at 1-4 km resolution, so the minimum spacing used here (2W = 0.5 km) would be too close together to resolve. However, we have observed adjacent cracks within the latter 2.5 km spatial scale. Similarly, in an H = 30 km ice shell at Europa, our fracture spacing ranges from 3 km (2W = 0.1H) to 15 km (2W = 0.5H), both of which are plausible given observed patterns on the surface. Maximum penetration depths for such fracture spacing range between 17% (0.1H = 3 km) and 28% (0.5H = 15 km), or ∼5.1 and 8.4 km, respectively, under 10 MPa tectonic stress. Comparing these depths to the single fracture case, a fracture reaching 43%-60% of a 5 km ice shell could occur at <2 MPa (5 times less); a 30 km ice shell would be fractured down to 17%-28% of its thickness under <5 MPa (2 times less). These estimates of stress are an order of magnitude less than the stress required in the multiple-fracture case.
We observe similar results for the range of ice shell thicknesses considered for Enceladus (Figures 3, 4). A 5 km ice shell would be fractured down to ∼90% of its thickness (4.5 km) if fractures were spaced 0.5H apart under 10 MPa tectonic stress. That same ice shell would likely be completely fractured under ∼2-3 MPa in the single fracture case. In the thickest ice shell considered, H = 90 km, if fractures are spaced 0.5H apart, they would penetrate to roughly 50% of the thickness (∼45 km) under 10 MPa; a single fracture could propagate to ∼45 km under ∼2 MPa. Figure 4 illustrates the higher stress required to fracture such fractured ice shells. Considering moderate stress values of 2-4 MPa, fractures spaced less than 0.5H in any Europa shell thickness cannot penetrate more than halfway through the ice shell. At Enceladus, under 2-4 MPa, no fracture spaced less than 0.5H apart from another fracture would penetrate farther than 80% of the way through the thinnest ice shell considered (e.g., the thickness at the south pole). In an ice shell thickness more relevant for Enceladus on the whole (60+ km), no fracture within that same spacing would propagate much farther than 50% of the thickness. There are no cases on either moon in which tidal stress (∼100 kPa) could fully fracture the ice shell in the presence of other fractures unless the ice was much less than a kilometer thick.

Isolated Basal Fractures
From the base of an ice shell, a single fracture filled with water is capable of penetrating further than a dry surface crack ( Figure 5). For example, at Europa, a single surface fracture propagating in a 5 km ice shell under tidal stress of 100 kPa would propagate to 125 m deep; if, instead, a fracture started from the base, a single water-filled fracture under the same opening stress would propagate to ∼1.3 km height. At Enceladus, a basal fracture under 100 kPa tidal stress in a 5 km ice shell could propagate to ∼4.7 km (94% of the total thickness). In fact, a single fracture at the base of Enceladus's ice shell, no matter the thickness (5-90 km), could propagate through the entire thickness (h/H > 95%) if the tectonic stress is at least 4-5 MPa. At Europa, a 5 km ice shell could be completely fractured (h/H > 95%) from the base under ∼2 MPa; a 30 km ice shell could be >95% fractured under 10 MPa. Comparing this with the single surface crack, it is clear that water-filled basal fractures dominate in terms of percentage fractured.

Multiple Adjacent Basal Fractures
Setting a basal fracture among adjacent fractures as in the surface case, Figure 6 similarly shows that a single isolated fracture will penetrate higher into an ice shell than in the multiplefracture case. Under the same 4-5 MPa load that would completely fracture (h/H > 95%) any Enceladus ice shell thickness (5-90 km) in the single basal fracture case, a 90 km ice shell would be fractured to between 60% and 73% of its thickness (54-65 km) depending on fracture spacing (0.1-0.5H), while a 5 km ice shell would be fractured to between 90% and 93% of its thickness (4.5-4.65 km). At Europa, while a single basal fracture would fully penetrate a 5 km ice shell under 2 MPa stress, that same stress would only drive a fracture 43%-68% through the same ice shell thickness, depending on spacing. While 10 MPa of stress could cause full rupture (h/H > 95%) of a 30 km Europa ice shell, the same stress would only drive a fracture  Figure 3). The error associated with these estimates from the polynomial fits (Section 2.2.1; not shown to reduce overlapping) is ±1% for all d/H values. Figure 5. Fractional penetration heights (h/H) of isolated water-filled basal fractures in Enceladus (left) and Europa (right) ice shells of varying thicknesses over a range of tectonic stresses. The base of the ice shell is assumed to be in hydrostatic equilibrium; thus, water at the base of the shell would rise through a crack to the hydraulic head (Section 2.2.2), applying pressure to crack walls that would act to open fractures. 52%-70% through the same ice shell thickness, again dependent on fracture spacing.

Shear Failure
In most cases investigated here, emplacement of adjacent fractures allowed neither surface nor basal crevasses to penetrate the entire ice thickness. Considering these results, we can turn to our discussion of stability against shear failure to determine scenarios in which full ice shell failure could be expected along optimally oriented faults above basal fractures. Figure 7 shows the minimum thickness of L S , the intact layer above a basal fracture, that is stable against shear failure under a range of tidal loads up to the expected peak 100 kPa (e.g., Hurford et al. 2007). For example, assuming a coefficient of friction of 0.3 in a 30 km Enceladus ice shell, the intact layer L S would have to be less than 6 km deep to allow for shear failure along optimally oriented faults under tidal stresses of 100 kPa. That is, a basal fracture would have had to propagate upward 24 km, to within 6 km of the surface, to allow the possibility for  . Minimum stable thickness of ice above a basal fracture in the ice shells of Enceladus (left) and Europa (right), across a range of tidal stress values (up to 100 kPa) oriented to cause shear slip. Solid lines and circles show the minimum thickness of ice that must be intact above a basal fracture in order for the unbroken part of the ice shell to be stable against shear failure under tidal stress, using a value of 0.3 for the coefficient of friction. Dashed lines of the same color show the range of results using different values for the coefficient of friction between 0.1 (lower bound) and 0.6 (higher bound). For a given ice shell thickness, if the intact layer is thicker than the stability boundary depth demarcated as solid lines, a basally fractured ice shell is stable against tidal shear failure. full ice shell rupture via shear failure. Similarly, in a 30 km Europa ice shell, a basal fracture would have to propagate upward to within ∼2.9 km of the surface to allow for shear failure and complete ice shell rupture. The effect of a higher acceleration due to gravity at Europa is at work here, enabling an overburden pressure that can accommodate the amplified tidal stress in the intact layer to a higher degree and allow for a closer approach by basal fractures to the surface without becoming unstable. The next question becomes, when do these scenarios allow for full ice shell failure through a combination of tensile and shear failure, i.e., how do these stability thickness boundaries (Figure 7) compare to the estimated basal fractures computed in Section 4.4 ( Figure 6)? Figure 8 shows the predicted maximum penetration height of basal fractures propagating at a distance of 0.1H (solid line) to 0.5H (dotted line of matching color) apart. These represent scenarios where fractures propagate at a distance between 1/10 and 1/2 of the ice shell thickness from each other. As shown in Figure 6, a spacing of 0.5H allowed for greater penetration than did an array spacing of 0.1H. Yellow, green, and blue circles in Figure 8 compare the minimum stable depth (from Figure 7) of L S at 150, 100, and 50 kPa, respectively. Comparison of these points to fracture heights shows the scenarios in which full ice shell failure could be achieved through a combination of tensile (Mode I) and shear failure.
For Enceladus (Figure 8, top), if tectonic stress at the base of the ice shell is near ∼1 MPa, basal fractures spaced half of the ice shell thickness (0.5H) apart could propagate to heights above the stability boundary depth in ice shells of thickness 15 km or less, and thus would be susceptible to shear failure under the tidal stress range. If the fractures were spaced more closely together, e.g., 0.1H (1/10 of the ice shell thickness), shear failure under peak tidal stress could be achieved if tectonic stress at the base of the ice shell reached ∼3 MPa. In a 30 km Enceladus ice shell, basal fractures can propagate up to the stability boundary L S if there is at least 1-4 MPa tectonic stress present to propagate them from the base in the 0.5H spacing case, or at least 3-8.5 MPa present at the base to propagate fractures in the 0.1H case. Fractures formed by stresses less than 1 MPa at the base in both cases would not reach the stability boundary depth; thus, the ice shell would be stable against shear failure under tidal stress. In a 60 km thick Enceladus ice shell, if tectonic stresses are present at the base of the shell that are at least 5 MPa or greater, basal fractures will propagate to within ∼12 km of the ice surface in the 0.5H spacing case, the stability boundary depth under 100 kPa tidal stress, allowing for the possibility of full ice shell rupture via shear failure. Otherwise, if fractures are spaced closer together in the 60 km ice shell (e.g., 0.1H), basal fractures do not propagate above the stability boundary depth under any stress considered (<10 MPa). In a 90 km Enceladus ice shell, we do not find any scenarios in which basal fractures are driven high enough in the ice shell to pass the stability boundary when considering fracture spacing between 0.1H and 0.5H, though this can be achieved in the single fracture model.
At Europa, as shown by the shallower stability boundary depths in Figure 7 compared to Enceladus, the enhanced overburden pressure enabled by the higher acceleration due to gravity requires higher fractures to be driven from the base of the ice shell (i.e., higher stresses required at the base). For a 5 km thick Europa ice shell, assuming tidal stresses of 100 kPa, the stability boundary depth is ∼2 km. Tensile basal fractures can reach this height if there is an applied tectonic stress at the base of at least 2.5 MPa in the 0.1H spacing case, or >1 MPa in the 0.5H spacing case. For a 10 km Europa ice shell, a tidal stress of 100 kPa results in a stability boundary depth of ∼2.8 km and basal fractures spaced at least 0.5H apart and can be propagated to this height under tectonic stresses of 5.5 MPa. If fractures are spaced 0.1H apart, there are no scenarios under which shear failure can connect basal fractures under tidal loading of 100 kPa. If 150 kPa additional loading were applied (above expected peak tidal stress; yellow circle), basal fractures propagated from the base under >7 MPa could potentially enable shear failure. The stability depth against tidal loading in a 20 km Europa ice shell is 4 km; in a 30 km ice shell this stability boundary is at ∼5 km. We do not find any scenarios in which basal fractures are driven high enough in either of these ice shell thicknesses to pass the stability boundary when considering fracture spacing between 0.1H and 0.5H, though this can be achieved in the single fracture model (see Figure A1).
Thus, we find that, assuming that large-enough basal fractures form in the subsurface under tectonic tensile stresses, further failure via tidal shear along optimally oriented faults above such basal fractures could lead to full ice shell rupture in shallow ice shell settings. For tectonic stress regimes between 1 and 3 MPa, for instance, only shells that are <5 km for Europa and less than 15 km for Enceladus are susceptible to complete ice shell rupture (by a combined tensile basal fracture and associated shear failure). These results were determined using an assumed coefficient of friction for ice of 0.3, a commonly applied value (Schulson & Duval 2009); however, this value in lab and modeling experiments has also ranged between 0.1 and 0.6 (Smith-Konter & Pappalardo 2008). The shaded regions in Figure 7 illustrate the effect of that range of friction values (μ f = 0.1 and 0.6). It can be noted in both cases that a higher frictional coefficient leads to a thinner layer being stable against shear, while a lower frictional coefficient demands a thicker section of intact ice above a basal fracture to remain stable against shear.

Discussion
Models of isolated fracture propagation in ice shells underestimate the magnitude of stress required to form surface-subsurface conduits. The effect of accounting for adjacent fractures is to lower the stress intensity as tensional stress is reduced in the slabs of ice separating them. Thus, we find that larger tectonic stresses are required to completely fracture through an ice shell. We argue that even subparallel fractures would behave in a similar manner, in that blocks of ice between fractures hold a reduced tensile stress. This suggests that as the ice shells became increasingly fractured throughout their evolution, it became difficult for new fractures to propagate deeper into the ice shell. Figure 9 shows the total fractional thickness of Enceladus (top) and Europa (bottom) ice shells expected to be fractured (basal and surface combined) using our tensile-only approach.
A particularly tempting set of fractures to investigate are the parallel tiger stripe fractures at Enceladus's south pole. The tiger stripes have a spacing W = 35 km; our results show that a 90 km ice shell with W = 35 km would not permit full ice shell fracture, nor would a 60 km ice shell. However, when 2W > H, the multiple-fracture consideration is invalid (e.g., Tada et al. 2000). Thus, using the thin-shell methodology (W = 35 km > H = 5-40 km; Equations (3) and (7)), under 2 MPa of stress, the tiger stripe fractures could have propagated through a 15 km ice shell, if they initiated from the base; in a 5 km ice shell, base-initiated tiger stripes could have reached the surface under ∼0.15 MPa of stress, on the order of tidal stresses there.
Our LEFM model predicts that closely spaced fractures are unlikely to penetrate the entire ice shell of Enceladus or Europa. However, fluid-filled basal fractures may increase the penetration of fractures within the deeper reaches of an ice shell. This possibility is supported by the prevalence of water-filled basal crevasses in terrestrial ice shelves (e.g., Luckman et al. 2012;McGrath et al. 2012a). Given that the base of the ice shell is likely near strengthless, after initial failure the basal fracture could be  Figure 6) over the range of ice shell thicknesses tested. Solid lines represent fracture spacing of 0.1H; dotted lines represent fracture spacing of 0.5H (same limits in Figure 6). Colored circles on each line show the scenarios in which the predicted fracture height reaches the minimum thickness for stability against shear failure. These are shown over a range of possible additional applied stress values (including max. tidal stress of 100 kPa (green) plus and minus 50 kPa (yellow, 150 kPa; blue, 50 kPa). Any thickness of ice over a fracture that is less than these values (i.e., any point on a line above and to the right of a circle) is susceptible to shear failure under tidal loading; any point on a line lower and to the left of a colored circle would not be susceptible to shear failure under tidal loading. For complete ice shell failure, both scenarios must be in play (high-enough tidal stress to cause shear, and a large-enough basal fracture to allow for it). expected to refreeze relatively quickly, save for additional changes to the stress field (e.g., Buffo et al. 2020). The expulsion of gas during this freezing process leaves weak bands that could also refracture under certain stress. Water pressure, for example, at least partially counteracts the closing effect of overburden pressure, further enabling the opening of a fracture.  Figures 3-6). Spread in color represents fracture spacing between 0.1H and 0.5H. While it is unlikely that surface and basal crevasses exactly line up in a real ice shell, it is likely that some such fractures exist close enough to each other in the horizontal dimension to be well aligned enough to enable low-dip angle faults to form between them. This plot provides an upper bound for the maximum fraction of the ice shell that can be fractured by a given tensile stress.
If a basal fracture were to penetrate a significant fraction of the ice shell thickness (Figure 8), shear failure above such a fracture may allow full ice shell rupture. This process has been invoked previously to define the stability of terrestrial ice cliffs . Support for this theory comes from the fact that in many geological materials, failure occurs along faults or slip lines when shear stress exceeds the yield stress (Turcotte & Schubert 2002). It remains to be determined whether shear failure occurs on Europa (Pappalardo et al. 2009), but we find that full ice shell rupture is plausible once a shell has been fractured via tensile stress to a certain fraction (Figures 7, 8). Thus, despite the difficulty of propagating a tensile fracture all the way downward between the surface and subsurface ocean, ice shells that are already riddled with basal fractures could plausibly be completely ruptured from the bottom up. This consideration requires the possibility that although some deep fractures could have resulted from large onset stresses (∼10 MPa) earlier in a planet's evolution (e.g., nonsynchronous rotation, true polar wander), those stresses are not necessarily the ones that caused, or are currently causing, full-shell fracture. It is similarly unlikely that tidal stress causes deep or full-shell fracture.
In Section 4.5, we discussed the penetration required by basal fractures (height h) in proximity to other fractures (spaced 0.1H to 0.5H apart) at either Enceladus or Europa to allow for shear failure of the ice shell via tidal stresses. This analysis showed that if initial tensile stresses were large enough in either shell, it is possible that tidal stresses have opened surfaceocean conduits along well-aligned surface and basal fractures. This again showcases that the processes that form large fractures are not necessarily the same ones that allow them to propagate. To allow for full penetration of the ice shell, optimally oriented faults are assumed to exist over basal fractures, and these fractures are assumed to be oriented in a direction that makes them susceptible to adequate tidal stress. The angle at which this tidal stress is applied has been shown to vary throughout an orbital period and may not activate fractures that are not optimally oriented to the direction of applied tidal stress. Thus, this is a lower bound on the possibilities for combined fracture modes to allow for full ice shell failure. That said, we found that ice shells between 15 and 30 km or less (depending on fracture spacing) at Enceladus could plausibly be fully ruptured under a combination of tectonic stresses at the ice shell base under 3 MPa (propagating basal tensile fractures) and stresses on the order of tidal stresses near the surface (peak 100 kPa) causing additional shear failure along faults above those basal features. At Europa, if basal tectonic stresses exist up to ∼3 MPa, ice shell thicknesses less than 5-10 km (depending on fracture spacing) could plausibly fully rupture if basal fractures can propagate up to within the stability boundary depth (Figure 8), and tidally induced shear stresses allowed for shear failure along optimally oriented faults above those fractures. Each of these values for both ice shells assumes that optimally oriented faults exist above the basal fractures. Of course, it is also possible that surface fractures as discussed in Section 4.2 could overlie these basal fractures and could theoretically serve to lessen that stable boundary depth L S . Thus, again, this analysis shows the lower bound for plausible full-shell rupture scenarios through a combination of tensile and shear faults activated or opened by diurnal tidal stress.
An additional implication of our study are fracture opening rates and source regions for possible plumes. On both bodies, fractures have been considered the origin of observed (Enceladus) and putative (Europa) plumes of gas and water (Porco et al. 2006;Hurford et al. 2007;Roth et al. 2014) under the control of tidal and normal stresses. Tidally varying periodicity has been confirmed for Enceladus, though there is a lag between tidal maximum and eruption rate (Hedman et al. 2013;Běhounková et al. 2015); as yet, no periodicity has been observed for Europa (e.g., Roth et al. 2014;Sparks et al. 2017). To accurately model the opening and closing of surface fractures, their setting among nearby fractures must be taken into account. Given that the south polar plumes are concentrated in 98 individual jets along the fractures, the fractures must be very narrow (Porco et al. 2006), i.e., tidal stresses are not capable of opening the tiger stripes to the extent previously thought. It is likely that once a fracture penetrates to the depth of a subsurface reservoir, enabling active venting, volatiles and gas filling the fractures lead to increased pressure on fracture walls (e.g., Kieffer et al. 2006). This increased pressure may lead to a similar effect of water pressure in basal fractures, in that they oppose the closing effect and enable further opening of fractures. Currently, it is particularly difficult to realistically constrain these effects on Europa, because the resolution of images obtained over the putative plume source regions is too low to determine fracture spacing there. That said, an area of future exploration would be to investigate conditions required for the formation and effects of plume eruption (from a high-pressure subsurface reservoir; e.g., Johnston & Montési 2017) on fracture opening rates. One possible driver of increased opening pressure on basal fracture walls could be the refreezing of the water in the fracture itself (Buffo et al. 2020).

Conclusions
Fractures that propagate in a field of fractures are less likely than isolated fractures to penetrate the full ice shell thickness, unless larger tensile stresses are present than previously reported. We found no scenarios in which surface or basal fractures, spaced closer together than half of the ice shell thickness, would propagate through the entire ice shell of any thickness tested under purely extensional stress. We have calculated envelopes of upper layer stability against shear failure above basal fractures, and our results advocate that the most likely means of enabling through-fracture of these ice shells is the formation of large basal fractures and subsequent shear failure along faults connecting those deep basal fractures to the surface in thinner ice shells (Europa, <10 km; Enceladus, <30 km). As such, it is likely that the processes that formed deep partial-shell fractures are not tidally induced; rather, they are related to the larger-scale stresses that existed earlier in these moons' histories. Our results firmly suggest that the fact that the shells are so fractured has significant impact on the feasibility of complete through-fracture to a deep subsurface ocean or reservoir. This conclusion serves to suggest that the highly fractured state of the ice shells at Europa and Enceladus should be considered in models of fracture propagation and in models of subsurface material escape. This work was funded over an elongated period of years by several grants and programs at different stages of the work, including partial support from the Michigan Space Grant Consortium, NASA grant NNX10AB216G, NSF CAREER-NSF-ANT 114085, and NSF grant EAGER-NSF-ARC 1064535. Some of the research was carried out and the manuscript was completed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, where it was supported by NASA's Cryosphere and Postdoctoral Fellowship Programs. Special thanks go to several advisors over the years, including Dr. Kate Craft for constructive support for the topic. Figure A1 shows results for the single fracture case for basal fractures in the ice shell of Enceladus and Europa. These are provided for comparison against the multiple fracture results in Figure 8 of the main text to illustrate the diminished capacity of the same applied stress to penetrate the ice shell when proximal cracks are taken into account. Figure A1. Potential shear failure scenarios in the ice shell of Enceladus (left) and Europa (right) for single basal fractures only (for comparison with Figure 8, main text, which shows the results in the multiple-fracture case). Lines show the minimum depth between the surface and the top of a basal fracture (i.e., distance between the surface and maximum height of fractures determined in Section 4.3, Figure 5, main text) over the range of ice shell thicknesses tested. Colored circles on each line show the scenarios in which the predicted fracture height reaches the minimum thickness for stability against shear failure. These are shown over a range of possible additional applied stress values, including max. tidal stress of 100 kPa (green) plus and minus 50 kPa (yellow, 150 kPa; blue, 50 kPa). Any thickness of ice over a fracture that is less than these values (i.e., any point on a line above and to the right of a circle) is susceptible to shear failure under tidal loading; any point on a line lower and to the left of a colored circle would not be susceptible to shear failure under tidal loading. For complete ice shell failure, both scenarios must be in play (high-enough tidal stress to cause shear, and a large-enough basal fracture to allow for it). The single isolated fracture case shown here, as opposed to the multiplefracture case in Figure 8, suggests that this would be a widespread mechanism to open fractures between the surface and the ocean.