Europa, one of Jupiter's icy moons, has features that look similar to mid‐ocean spreading ridges, and recent observations suggest that Europa may also have something similar to subduction zones. This implies that Europa may be the only planetary body other than Earth to have plate tectonics. We construct a computer simulation to track the temperature and compaction of ice as it sinks into warmer underlying ice. When ice warms up, it becomes less dense and easier to compact. The salt content of Europa's plates is also important because salt is denser than ice. The temperature, salt content, and how compacted the ice is determine the density and whether the ice can sink. We find that if the salt content of Europa's ice shell varies by few percent from place to place, then the plates can sink. This supports the idea that something like plate tectonics is occurring on Europa and may tell us about the composition of Europa's ice shell. Our work also implies that the plates will sink all the way to Europa's subsurface ocean. This is important because material from the surface of Europa could act as food for life that may exist in Europa's ocean.

Motivated by recent evidence for subduction in Europa's ice shell, we explore the geophysical feasibility of this process. Here we construct a simple model to track the evolution of porosity and temperature within a slab that is forced to subduct. We also vary the initial salt content in Europa's ice shell and determine the buoyancy of our simulated subducting slab. We find that porosity and salt content play a dominant role in determining whether the slab is nonbuoyant and subduction in Europa's ice shell is actually possible. Generally, we find that initially low porosities and high salt contents within the conductive lid are more conducive to subduction. If salt contents are laterally homogenous, and Europa has a reasonable surface porosity of ϕ 0 = 0.1, the conductive portion of Europa's shell must have salt contents exceeding ~22% for subduction to occur. However, if salt contents are laterally heterogeneous, with salt contents varying by a few percent, subduction may occur for a surface porosity of ϕ 0 = 0.1 and overall salt contents of ~5%. Thus, we argue that under plausible conditions, subduction in Europa's ice shell is possible. Moreover, assuming that subduction is actively occurring or has occurred in Europa's recent past provides important constraints on the structure and composition of the ice shell.

1 Introduction The surface of Europa, one of Jupiter's icy satellites, is marked by a myriad of extensional tectonic features (e.g., Greenberg et al., 1998; Hoppa et al., 1999). The dearth of craters on Europa implies a young surface age of 40–90 Myr and that Europa is actively resurfaced (Bierhaus, Zahnle, & Chapman, 2009). The ubiquity of extensional tectonics, including dilational bands similar to mid‐ocean spreading ridges (Prockter et al., 2002; Sullivan et al., 1998), and the lack of corresponding compressional tectonics have led to the hypothesis that resurfacing on Europa is achieved by subduction or subsumption of portions of the ice shell (Kattenhorn & Prockter, 2014). Detailed plate reconstructions act as strong support for this hypothesis (Kattenhorn & Prockter, 2014). Similar plate reconstructions focused on other regions of Europa also find evidence of subduction (Collins et al., 2016; Perkins, Patterson, & Prockter, 2017). As a descending plate warms up, it may be integrated into warmer underlying ice, or subsumed. If, however, the plate remains dense and continues to descend, subduction is a more apt term. We will use the term subduction in this paper and later argue why this term may be more appropriate. It is unclear if subduction is initiated and driven by convection (Kattenhorn & Prockter, 2014) or if other driving forces such as nonsynchronous rotation (Geissler et al., 1998; Wahr et al., 2009) or true polar wander (Matsuyama & Nimmo, 2008; Schenk, Matsuyama, & Nimmo, 2008) aid in driving plate motions and subduction. When including plasticity and yield strengths less than 100 kPa, simulations of convection in Europa's ice shell produce conductive lids and surfaces that are continually overturning or have episodic catastrophic overturns (Showman & Han, 2005). O'Neill and Nimmo (2010) similarly advocate for a relatively low yield strength to explain episodic overturn and high heat flows in the south polar terrain of Saturn's moon Enceladus. Despite the overturn and recycling of the surface occurring in some Europan convection models (Showman & Han, 2005), we identify major questions about the geophysical feasibility of subduction in Europa's ice shell. The density difference between solid ice and a liquid ocean makes it essentially impossible to subduct a purely conductive ice shell into an underlying ocean. However, a colder conductive ice shell could subduct into a warmer underlying ice layer, since colder ice is denser and would therefore be negatively buoyant relative to surrounding warm ice. Thus, as did Kattenhorn and Prockter (2014), we assume that the conductive portion of the ice shell is subducting into a warm convecting portion of the ice shell (e.g., Nimmo & Manga, 2009) (Figure 1). For a conductive lid with a surface temperature of 100 K, a basal temperature of 260 K, and a thermal conductivity proportional to 1/T, the typical thermal anomaly between the subducting conductive lid and underlying warm convecting ice will be ΔT≈93 K (Nimmo, Pappalardo, & Giese, 2003). The density anomaly associated with this temperature difference is given by Δρ = − αρΔT, where α = 1.56 × 10−4 (T/250 K) K−1 is the temperature‐dependent volumetric coefficient of thermal expansion and ρ = 920 kg/m3 is a reference density typical for ice (Kirk & Stevenson, 1987). The average density anomaly is thus Δρ≈ 9 kg/m3, where the subducting conductive lid is denser than the warmer underlying ice. This simple analysis does not, however, consider porosity in the conductive shell. If the conductive portion of the ice shell is on average more than 1% more porous than the underlying warm convecting ice, the “subducting” shell will be buoyant (the density anomaly associated with porosity differences, Δϕ, is given by Δρ = − ρΔϕ). It is difficult to imagine that subduction could occur if the shell were buoyant. Figure 1 Open in figure viewer PowerPoint h cond , and assumed viscosities at the base of the conductive portion of the ice shell η b (solid curves η b = 1015 Pa s; dashed curve η b = 1013 Pa s). Range of (a) possible temperature and (b) normalized porosity within the slab at the moment subduction begins (time = 0). Initial temperature (Figure 1 a) and normalized porosity (Figure 1 b) plotted as a function of depth for different conductive shell thicknesses,, and assumed viscosities at the base of the conductive portion of the ice shell(solid curves10Pa s; dashed curve10Pa s). The porosity at Europa's optical surface may be higher than 95% (Hendrix, Domingue, & King, 2005). The 70 cm radar scattering of Ganymede and Callisto, the other two icy Galilean satellites, imply near surface porosities exceeding 40% (Black, Campbell, & Ostro, 2001). Modeling of 3.5 cm, 13 cm, and 60 cm wavelength radar reflectivities suggests that Europa has a surface layer with appreciable porosity that is meters thick (Black, Campbell, & Nicholson, 2001). Although Europa's near surface has significant porosity, little is known about porosity in Europa's ice shell at depths exceeding a few meters. Recent results of NASA's Gravity Recovery and Interior Laboratory (GRAIL) mission (Zuber et al., 2013) suggest that the Moon's crust is 12% porous to depths of at least a few kilometers (Besserer, Nimmo, & Wieczorek, 2014; Wieczorek et al., 2013). Although porosity on an icy satellite may be quite different than porosity on the Moon, pore space is likely created by fracturing and dilatancy, increase in rock volume with deformation (e.g., Jaeger, Cook, & Zimmerman, 2009), as the ice shell is tectonized. Thus, the extensional tectonic features covering the surface of Europa (e.g., Greenberg et al., 1998; Hoppa et al., 1999) suggest that the conductive portion of Europa's ice shell may have appreciable porosity. Of course, the porosity and temperature of the “subducting” ice will evolve with time. As the shell subducts, pore space in the ice will eventually collapse and diffusion of heat will erase the thermal anomaly. Because the rate of pore space crushing depends on material viscosity, and thus temperature, the evolution of porosity will be linked to the thermal evolution of the subducting slab. This coupled evolution of porosity and temperature in a subducting ice shell is the focus of this work. In section 2 we describe our numerical model and how we track the evolution of porosity and temperature of a subducting slab. In section 3.1 we describe our fiducial model and the effect that initial porosity has on the buoyancy of the subducting slab. In section 3.2 we study the role of the temperature structure within Europa's ice shell. In section 3.3 we discuss the effects of bending, subduction dynamics, and uncertainties in the rheology of ice. In section 3.4 we explore the consequence of initial salt content, including lateral variations and the conditions needed to make subduction possible. Finally, in section 4 we discuss implications of our work for the structure and composition of Europa's ice shell.

2 Model Description We use a simple one‐dimensional model to track the evolution of temperature and porosity within a slab as it subducts through a reference ice shell, consisting of a conductive lid and an underlying warm convective layer (Figure 2). Figure 1 shows the range of initial (time = 0) porosity and temperature profiles of the slab used in our models, and Figure 2 presents a schematic of the model. Our simulations track the evolution of an individual parcel or slice of the slab. The parcel is broken into 200 equal‐sized cells. At each time step, the depth of each cell is evolved assuming a constant slab velocity and a subduction geometry, as shown in Figure 2. As the depth of the parcel increases, the top temperature boundary condition changes and overburden pressure increases. At each time step, temperature is evolved using Fourier's law of heat conduction and porosity is evolved by tracking removal of pore space by viscous flow. We discuss the model in more detail below. Figure 2 Open in figure viewer PowerPoint Geometry of a simple one‐dimensional model of a subducting slab. The modeled parcel of material is shown at three times, time = 0, t bend, and t n . Note that the schematic only shows 5 cells, while our actual model has 200 cells. Additionally, for a time step Δt = 1 year and spreading rate v spread = 40 mm/yr, the cells have a horizontal extent Δx ≈ v spread Δt = 2.8 cm while for h cond = 6 km the vertical dimension Δy = h cond /n cell = 30 m, where n cell is the number of cells. The dark blue material at the top left represents the dilational spreading bands. The schematic shows a slab bent through the dip angle θ = 45° along a radius R bend ≈ 2 h cond in a time t bend = R bend θ/v spread . After time t bend , slab continues to subduct at a constant dip angle and constant rate, v spread . τ given by (1) First, we consider whether a one‐dimensional model is a reasonable representation of a subducting slab. Thermal density differences will persist for a timescalegiven by (Turcotte & Schubert, 2002), where the thermal diffusivity κ = 10−6 m2 s−1 and h cond is the typical length scale of the thermal anomaly (Klinger, 1980; Nimmo et al., 2003). For h cond = 6 km, the thickness of the conducting portion of the ice shell in our fiducial model, τ ≈ 1 Myr. The distance the subducting slab will travel during this time is given by l = v spread τ, where v spread is velocity of the subducting slab. A comparison of the thickness of the conducting portion of the shell, h cond (i.e., the thickness of the subducting slab in our model), and this distance, l, provides an estimate of the importance of thermal diffusion along the direction of slab motion. For example, if l ≈ h cond , heat diffusion from the top to the bottom of the slab will be comparable to the lateral diffusion of heat along the slab. If, however, l ≫ h cond , heat will primarily diffuse through the top of the slab and a one‐dimensional model is reasonable. To compare l and h cond , we first need an estimate of the velocity of the subducting slab, v spread . Following Kattenhorn and Prockter (2014), we assume that the slab velocity is equal to the estimated spreading rate at Europa's dilational bands. Based on a mid‐ocean spreading ridge model and a relationship between fault block spacing and v spread , Stempel, Barr, and Pappalardo (2005) argue that v spread could be up to 40 mm/yr. Kattenhorn and Prockter (2014) suggest that this rate is consistent with subduction on Europa being a dominant resurfacing mechanism. Convection models with overturning lids predict even higher surface velocities of 100–1000 mm/yr (Showman & Han, 2005). At a spreading rate of 40 mm/yr, v spread τ≈46 km, substantially larger than the thickness of the subducting slab. Thus, a one‐dimensional conduction model will provide reasonable estimates of temperature evolution within a subducting slab. T s = 100 K and temperature at the base of the conducting portion of the slab of T b = 260 K. The thermal conductivity of ice is given by Petrenko and Whitworth ( 1999 (2) We assume an initial surface temperature ofK and temperature at the base of the conducting portion of the slab ofK. The thermal conductivity of ice is given by Petrenko and Whitworth () as T, the steady state temperature of the conducting portion of the shell is given by (3) z is depth and h cond is the assumed thickness of the conductive portion of the shell ( h cond = 6 km, in our fiducial model) (Nimmo et al., 2003 h cond the temperature is assumed to be constant at 260 K (Figure 2003 2002 For a temperature‐dependent thermal conductivity, proportional to 1/, the steady state temperature of the conducting portion of the shell is given bywhereis depth andis the assumed thickness of the conductive portion of the shell (km, in our fiducial model) (Nimmo et al.,). Belowthe temperature is assumed to be constant at 260 K (Figure 1 ) approximating the behavior of nearly isothermal convecting ice (Nimmo et al.,). In an actual ice shell the temperature would increase in a thin boundary layer near the interface of the ice shell and ocean (e.g., Hussmann, Spohn, & Wieczerkowski,). ϕ 0 , ice shell that is allowed to evolve for t por = 65 Myr, in the middle of estimates of 40–90 Myr for the age of Europa's surface (Bierhaus et al., 2009 ϕ 0 depends on depth, so little is known about porosity at depth in Europa's ice shell that we argue that the prudent choice is a constant ϕ 0 . This can be thought of as a depth‐averaged initial porosity, which is then allowed to evolve through time t por . This reference porosity as a function of depth is described exactly by (4) P = ρgz , g = 1.31 m/s2 and η(z) is the viscosity as a function of depth. The rheology of ice is often written as an effective Newtonian viscosity (e.g., Nimmo et al., 2003 (5) Q = 50 kJ/mol (Goldsby & Kohlstedt, 2001 η b = 10 15 Pa s (for our fiducial model) (Nimmo et al., 2003 R is the ideal gas constant, and T(z) is the temperature as a function of depth (Figure T can be substituted for T(z) in equation η(T) . Our results are rather insensitive to the choice of t por since viscosity η(z) is strongly dependent on temperature and depth. Thus, we do not expect a large change from the results presented here even if the subducting material is much older than 65 Myr or if fracturing and dilatancy associated with extensional tectonics (e.g., Greenberg et al., 1998 1999 t por . We assume the starting porosity of the ice shell as given in Figure 1 . This porosity structure represents an initially constant porosity,, ice shell that is allowed to evolve forMyr, in the middle of estimates of 40–90 Myr for the age of Europa's surface (Bierhaus et al.,). Although one could assume thatdepends on depth, so little is known about porosity at depth in Europa's ice shell that we argue that the prudent choice is a constant. This can be thought of as a depth‐averaged initial porosity, which is then allowed to evolve through time. This reference porosity as a function of depth is described exactly bywhere the pressurem/sandis the viscosity as a function of depth. The rheology of ice is often written as an effective Newtonian viscosity (e.g., Nimmo et al.,) given bywhere the activation energy iskJ/mol (Goldsby & Kohlstedt,), the reference viscosity of the base of the shell isPa s (for our fiducial model) (Nimmo et al.,),is the ideal gas constant, andis the temperature as a function of depth (Figure 1 ). Note thatcan be substituted forin equation 5 to give the temperature‐dependent viscosity. Our results are rather insensitive to the choice ofsince viscosityis strongly dependent on temperature and depth. Thus, we do not expect a large change from the results presented here even if the subducting material is much older than 65 Myr or if fracturing and dilatancy associated with extensional tectonics (e.g., Greenberg et al.,; Hoppa et al.,) create porosity on a time scale much shorter than We refer to the infinitesimal slice of the subducting slab that our model represents as a parcel of the slab. At time zero in our model, the center of the slab moves at a constant rate v spread along a radius of curvature R bend = 2 h cond (in our fiducial model) so that in a time t bend = R bend θ/v spread , the parcel has rotated through the dip angle θ (Figure 2). The depth of each cell is evolved according to this geometry (Figure 3). After t bend , the depth of each of the 200 cells is evolved according to the spreading rate and dip angle of the slab, assumed to be θ = 45° in our fiducial model. The dynamics of a bending plate and dip angle are essentially unconstrained on Europa. In section 3.3, we explore the effect of assumed bending dynamics by varying spreading rate, v spread , R bend , and θ. Although our model is quite simple, it serves as a starting point to explore the geophysical feasibility of subduction in Europa's ice shell. Figure 3 Open in figure viewer PowerPoint h cond = 6 km, v spread = 40 mm/yr, η b = 1015 Pa s, and an initial temperature and porosity structure as describe in Figure ϕ 0 = 0.1 . Evolution of (a) temperature, (b) porosity, and (c) depth in a subducting slab. Material in the slab is colored according to temperature (Figure 3 a) and porosity (Figure 3 b) as described by the color scale bars. These results are from our fiducial modelkm,mm/yr,Pa s, and an initial temperature and porosity structure as describe in Figure 1 , where (6) z top is the depth of the top of the parcel, T ref (z top ) is the temperature in the reference shell (Figure z top , and T 1 is the current temperature of the topmost cell. In practice, T 1 ≈ T ref and the boundary temperature is very close to T ref . Because the temperature of the reference ice shell T ref is assumed to be constant with time at a given depth, this boundary condition assumes that the cold subducting slab does not alter the temperature of the surrounding material and tends to overestimate heating of the subducting slab. One way to quantify that the effect of assuming T ref is constant with time at a given depth is to assume that the top cell of the parcel does not heat up at all. In this case, the top temperature boundary condition as a function of depth would be given by (7) The temperature boundary condition at the bottom of the slab is assumed to be a constant 260 K. The temperature at the contact of a warmer and cooler object of the same composition will be the average of their respective temperatures. Thus, the top temperature boundary condition of our model,whereis the depth of the top of the parcel,is the temperature in the reference shell (Figure 1 ) at a depth, andis the current temperature of the topmost cell. In practice,and the boundary temperature is very close to. Because the temperature of the reference ice shellis assumed to be constant with time at a given depth, this boundary condition assumes that the cold subducting slab does not alter the temperature of the surrounding material and tends to overestimate heating of the subducting slab. One way to quantify that the effect of assumingis constant with time at a given depth is to assume that the top cell of the parcel does not heat up at all. In this case, the top temperature boundary condition as a function of depth would be given by If the temperature of the reference ice shell was permitted to evolve, we would expect the top boundary temperature to be intermediate of the results given by equations 6 and 7 (note that equation 7 is only used for comparison; equation 6 is the boundary condition we actually use in our model). The extreme and clearly unphysical assumption for equation 7 would result in at most a 30% difference from the temperature boundary condition we used. Using a one‐dimensional model neglects the diffusion of heat along the direction of slab motion. In this way, the one‐dimensional model tends to underestimate heating of the slab. This may counteract the overestimate of heating caused by assuming that the temperature of the reference ice shell is unchanging. Thus, our assumed top boundary condition and use of a one‐dimensional model are reasonable for our order of magnitude estimates of slab buoyancy. (8) q is the heat flow per unit area, k is the temperature dependent thermal conductivity (equation y is the distance from the top of the parcel (Figure T(y) is the temperature as a function of y as calculated in the previous time step. The temperature is then evolved using the heat flow calculated from equation c = 1925 (T/250 K) J/kg/K (Kirk & Stevenson, 1987 2002 The evolution of our model through time is achieved through a simple explicit finite difference scheme. At each time step we calculate the heat flow in our model using Fourier's law of heat conductionwhereis the heat flow per unit area,is the temperature dependent thermal conductivity (equation 2 ),is the distance from the top of the parcel (Figure 2 ), andis the temperature as a function ofas calculated in the previous time step. The temperature is then evolved using the heat flow calculated from equation 8 and assuming a temperature‐dependent specific heat capacity for ice ofJ/kg/K (Kirk & Stevenson,). With 200 cells and temperature‐dependent thermal conductivity, we find that a time step of 1 year produces models that reliably converge. The temperature evolution model was validated against the analytical solution for instantaneous heating of a half‐space (Turcotte & Schubert,). ϕ , is removed by viscous flow and evolves with time at a rate depending on confining pressure and temperature‐dependent viscosity (9) 1985 η is the temperature dependent viscosity described by equation P = ρgz , and z is the current depth of the cell. The porosity ϕ and viscosity η in equation 2003 During each time step the porosity in each cell,, is removed by viscous flow and evolves with time at a rate depending on confining pressure and temperature‐dependent viscositygiven by Fowler (), whereis the temperature dependent viscosity described by equation 5 , andis the current depth of the cell. The porosityand viscosityin equation 9 are the porosity and viscosity from the previous time step. The porosity evolution model was validated using the analytical solution for pore space crushing in an initial constant porosity slab (i.e., equations 4 and 5 ) (Nimmo et al.,). After the porosity and temperature are updated, the depth of each cell in the parcel is updated according to the geometry shown in Figure 2. Increasing the depth results in an updated overburden pressure P = ρgz and an updated top temperature boundary condition defined by equation 6 and the current depth of the parcel. (10) ΔT , Δϕ , and Δf salt are the difference in temperature, porosity, and salt content, respectively, between the cell of interest and material in the reference ice shell at the same depth as the cell of interest (note that ΔT , Δϕ , and Δf salt are not the changes in temperature, porosity, and salt content with time). The temperature‐dependent volumetric coefficient of thermal expansion is given by α = 1.56 × 10−4 (T/250 K) K−1 and ρ = 920 kg/m3 is a reference density typical for ice (Kirk & Stevenson, 1987 ρ salt = 1444 kg/m3 appropriate for natron, a potentially dominant component of surface salts based on spectral unmixing analyses of McCord et al. ( 1999 2010 2016 Δf salt . The average density difference of the parcel, Δρ av , is then obtained by summing the density difference of each cell and dividing by the total number of cells. When Δρ av > 0 the parcel of the slab is denser than reference material. For reference, if the top of the parcel were instantaneously brought to a depth of h cond (meaning there is no change in temperature or porosity), Δρ av would be ~10 kg/m3, assuming ϕ 0 = 0 , Δf salt = 0 , and h cond = 6 km. For the same instantaneous comparison and ϕ 0 = 0.1 , Δρ av ≈ −41 kg/m3 resulting in a parcel that is buoyant. In addition to temperature and porosity, we also consider how salt content may affect the buoyancy of the parcel. To determine if the conducting parcel of the slab is buoyant, we compare the density of the parcel to the density of the reference ice shell at the same depth. We calculate the density difference of each cell using the following equationwhere, andare the difference in temperature, porosity, and salt content, respectively, between the cell of interest and material in the reference ice shell at the same depth as the cell of interest (note that, andare not the changes in temperature, porosity, and salt content with time). The temperature‐dependent volumetric coefficient of thermal expansion is given byandkg/mis a reference density typical for ice (Kirk & Stevenson,). We assume a salt density1444 kg/mappropriate for natron, a potentially dominant component of surface salts based on spectral unmixing analyses of McCord et al. (). Other spectral unmixing analyses suggest that magnesian sulfates dominate the surface (Shirley et al.,; Shirley, Jamieson, & Dalton,). As magnesian sulfates are more dense than natron (up to approximately twice as dense as natron), assuming that natron is the dominant salt component may produce conservative estimates for the density anomaly produced by a given difference in salt content,. The average density difference of the parcel,, is then obtained by summing the density difference of each cell and dividing by the total number of cells. Whenthe parcel of the slab is denser than reference material. For reference, if the top of the parcel were instantaneously brought to a depth of(meaning there is no change in temperature or porosity),would be ~10 kg/m, assuming, andkm. For the same instantaneous comparison and−41 kg/mresulting in a parcel that is buoyant.

3 Results 3.1 Fiducial Model and the Role of Porosity Figure 3 shows the results of our fiducial model with h cond = 6 km, v spread = 40 mm/yr, and an initial temperature and porosity structure as described in Figure 1, where ϕ 0 = 0.1. In this model, material in the slab remains colder than the underlying convecting ice into which it is being subducted for ~500 kyr. This timescale is shorter than implied by equation 1, because here we use a temperature‐dependent thermal conductivity and heat capacity. The porosity in the slab is completely crushed out by ~260 kyr. At v spread = 40 mm/yr, the parcel bends through the dip angle in a time, t bend ≈ 236 kyr. The slab then descends at a constant angle, and 100 kyr corresponds to increasing in depth by ~2.8 km. The average density anomaly of the parcel as a function of time (Figure 4) is crucial to understanding how the evolution of temperature and porosity affects the buoyancy of the slab. As we show many plots of average density anomaly versus time, it is prudent to discuss the important aspects of this plot. When the average density anomaly is positive, the modeled parcel of material is nonbuoyant, and we argue that subduction is feasible. When the average density anomaly is negative, the parcel of material is buoyant and subduction is unlikely. In our models, the parcel is forced to subduct with a depth‐time history as shown in Figure 1c (when R bend and v spread are unchanged from their values in the fiducial model). Many of our models result in buoyant parcels of material. By exploring the effect of various model parameters, however, we can determine the conditions that produce slabs that are less buoyant and more conducive to subduction. Using the information gained from sensitivity tests performed in other sections and including the possible role of salt content, in section 3.4 we demonstrate that plausible assumptions can produce conditions where subduction is feasible (i.e., nonbuoyant slabs). Figure 4 Open in figure viewer PowerPoint Effect of initial porosities, ϕ 0 , on the time evolution of the average density anomaly in the modeled parcel of subducting material. Different curves represent different initial porosities, ϕ 0 , in both the subducting slab and reference shell, as indicated by the legend. All other parameters are unchanged from our fiducial model. It is instructive to first focus on the density anomaly and time evolution of only the thermal anomaly, ΔT (i.e., ϕ 0 = 0; Figure 4, blue curve). As the parcel of material subducts, the temperature difference between each cell in the parcel and material at the corresponding depth in the reference model increases with time, as does the average density anomaly. The maximum density anomaly is achieved at ~240 kyr. As the parcel continues warming, this density anomaly decreases (Figure 4) until it ultimately disappears at ~500 kyr. For an initial porosity ϕ 0 = 0.1, the evolution of porosity dominates the evolution of the average density anomaly with time (Figure 4, black curve). In this case, the porosity difference, Δϕ, between the subducting slab and the reference material at the same depth increases with time resulting in the parcel of material being buoyant. The initial porosity is rapidly reduced at a depth of 3 km (Figures 1 and 3). Thus, when the modeled parcel of the slab reaches a depth of ~3 km, the difference in porosity between the slab and the reference material reaches a maximum (the top of the slab reaches a depth of 3 km slightly after 190 kyr when v spread = 40 mm/yr). After 200 kyr, crushing of pore space begins to dominate and pore space is completely removed by ~260 kyr (Figure 3). After this time, the density anomaly is well described by the thermal evolution and all curves converge (Figure 4). Between ϕ 0 = 0 and 0.1, the evolution is roughly a linear combination of these two end‐members. Motivated by Europa's high surface porosity (Black, Campbell, & Nicholson (2001); Black, Campbell, & Ostro (2001); Hendrix et al., 2005) and the results of NASA's GRAIL mission (Besserer et al., 2014; Wieczorek et al., 2013), we consider ϕ 0 = 0.1 as our fiducial model. The fiducial model is used as a reference to help understand the effect of changing various model parameters. Models with lower porosity tend to be less buoyant and more favorable for subduction than those with higher porosity. For porosity, less than ϕ 0 = 0.01, the parcel is nonbuoyant at all times and subduction is feasible (Figure 4). 3.2 Role of Temperature Structure The thickness of the conducting portion of Europa's ice shell is relatively uncertain (e.g., Nimmo & Manga, 2009). Figure 5 shows how our model results depend on the assumed thickness of the conducting part of Europa's ice shell, h cond . To isolate the thermal effects, we set the bending radius R bend = 12 km and have material at a depth of 3 km (i.e., the bottom of 3 km thick slab; middle of 6 km thick slab, and a third of the way down into the 9 km thick slab) move at a constant rate v spread . This is the same as the fiducial model and ensures that regardless of h cond , the slab bends toward the ultimate dip angle in the same amount of time, t bend . This also ensures that the evolution of the depth of the top of the parcel and the top temperature boundary condition with time remains the same as we vary h cond . We discuss the role of bending dynamics separately in section 3.3. Figure 5 Open in figure viewer PowerPoint Effect of the thickness of the conductive portion of the shell, h cond , on the time evolution of the average density anomaly in the modeled parcel of subducting material. Different curves represent different conductive layer thicknesses, h cond , in both the subducting slab and reference shell, as indicated by the legend. All other parameters are unchanged from our fiducial model. Because the surface temperature (T s = 100 K) and basal temperature (T B = 260 K) do not depend on h cond , a thicker ice shell corresponds to a lower heat flux and is colder at a given depth than a thinner shell (Figure 1a). The depth at which initial porosity becomes negligible also depends on h cond through equations 4 and 5 (Figure 1b). The mass anomaly created by an increasing difference in porosity with time is nearly the same in these different models. This mass anomaly, however, is averaged over a smaller thickness for models with smaller h cond resulting in an average density anomaly that becomes larger in magnitude at earlier times for smaller h cond . The thickness of the porous layer is nearly proportional to h cond (Figure 1b). Consequently, the minimum average density anomaly occurs later as h cond increases. This is because Δϕ is maximized at greater depth for larger h cond since the initial porosity extends to greater depths. We also note that the depths of the minimum average density anomalies are quite insensitive to h cond (Figure 5). These models clearly demonstrate that it takes longer for pore space to be crushed out of the thicker colder subducting slabs, as we might expect. We also see that the thermal anomaly disappears much faster in the thinner, warmer conducting slabs. For h cond = 3 km, the thermal anomaly disappears by ~250 kyr (Figure 5), while for h cond = 6 km the thermal anomaly persists to ~500 kyr (Figure 2). Although h cond affects the time evolution of slab buoyancy, the lack of a strong effect on the magnitude of buoyancy of the slab seems to suggest the feasibility of subduction as defined in this work is insensitive to the thermal structure of the shell. 3.3 Role of Ice Rheology and Dynamics of Subduction Following Nimmo et al. (2003), we also consider a range of basal shell viscosities, η b , from 1013 to 1015 Pa s to account for uncertainties in the rheology of ice. At the same temperature, pore space crushes out 100 times faster when η b = 1013 Pa s as compared to η b = 1015 Pa s (equation 9). Setting η b = 1013 Pa s leads to pore space completely crushing out by 220 kyr, whereas it takes 260 kyr when η b = 1015 Pa s (Figure 6). Lowering η b also reduces the depth to which the initial porosity extends (Figure 1), resulting in a lower magnitude for the minimum average density anomaly. These effects together make the slab less buoyant. Thus, models with lower η b are more conducive to subduction. Figure 6 Open in figure viewer PowerPoint Effect of the basal viscosity, η b , on the time evolution of the average density anomaly in the modeled parcel of subducting material. Different curves represent different basal viscosities, η b , in both the subducting slab and reference shell as indicated by the legend. All other parameters are unchanged from our fiducial model. The spreading rate, v spread , on Europa is another poorly constrained parameter (Showman & Han, 2005; Stempel et al., 2005). Increasing v spread forces the slab to subduct at a faster rate, meaning the parcel will reach greater depths at earlier times. Thus, a higher v spread results in larger overburden pressures and higher top temperature boundary condition (equation 6) at any given time. We explore the effect of v spread on the feasibility of subduction by varying v spread from 10 to 100 mm/yr (Figure 7). Note that Figure 7 has time on a logarithmic scale to accommodate the differences in timescale associated with this large range in v spread . Although the initial porosity structure is the same in these models (Figure 1b), we see marked differences in the initial slope of the average density anomaly as a function of time. A steeper slope at higher spreading rates results in a greater mismatch in porosity Δϕ at earlier times. For a higher spreading rate, the thermal anomaly disappears faster and pore space is completely crushed out at earlier times compared to models with a lower spreading rate. This faster pore space removal is explained by the top temperature boundary condition, which depends on depth, changing faster with the higher spreading rate. We note that this effect is not captured in a simple application of equation 1 for the lifetime of the thermal anomaly. As v spread decreases from 40 mm/yr to 10 mm/yr, the magnitude of the minimum average density anomaly decreases slightly. Thus, relatively lower spreading rates produce slabs that are less buoyant and slightly more conducive to subduction. Figure 7 Open in figure viewer PowerPoint Effect of the spreading rate, v spread , on the time evolution of the average density anomaly in the modeled parcel of subducting material. All other parameters are unchanged from our fiducial model. To further test the sensitivity of our models to the dynamics of subduction, we vary the radius of curvature the slab bends through, R bend . As R bend decreases, cells reach a greater depth at earlier times. Increasing R bend has a similar effect as reducing v spread , resulting in a lower magnitude minimum average density anomaly, which is reached at later times. Both the differences in timescale and the minimum average density anomaly are more sensitive to variation of v spread considered than to changes in R bend (Figures 7 and 8). Figure 8 Open in figure viewer PowerPoint Effect of the radius of curvature through which the slab bends, R bend , on the time evolution of the average density anomaly in the modeled parcel of subducting material. All other parameters are unchanged from our fiducial model. Finally, we vary the ultimate dip angle of the slab, θ (Figure 9). The ultimate dip angle of 30, 45, and 60° is reached at a time of 157, 236, and 314 kyr, respectively. The evolution of depth of the parcel is identical for θ= 45 and 60° up until 236 kyr. Because the pore space has been removed from the slab by 236 kyr, the evolution of parcel density is very similar for θ= 45 and 60° (Figure 9). Figure 3c shows that the depth of the bottom of the parcel increases faster once the ultimate dip angle has been reached and bending has ceased. This explains why the evolution of the average density anomaly for θ= 30° looks similar to θ= 45° with a slight shift to earlier times. Ultimately, Figure 9 shows that our results are insensitive to changes in θ. Figure 9 Open in figure viewer PowerPoint Effect of ultimate dip angle of the slab, θ , on the time evolution of the average density anomaly in the modeled parcel of subducting material. All other parameters are unchanged from our fiducial model. 3.4 Role of Salt Content Next, we consider the role of salts within Europa's ice shell. Radiative transfer modeling of visible and near‐infrared spectral data from the NIMS instrument aboard Galileo indicates that salt content on the surface of Europa is heterogeneous and may reach nearly 100% at some locations (McCord et al., 1999). Temperature, grain size, porosity, surface asperities or surface roughness, low albedo components, and trace contaminants, however, can greatly affect a model's salt content estimates (Dalton et al., 2005; Hapke, 1981, 1984; McCord et al., 1999; Milliken & Mustard, 2007a, 2007b). Nimmo et al. (2003) argue that warm ice will be relatively pure because salts will be incorporated in eutectic liquids that drain out of the warm, convecting ice on short timescales. Thus, we might expect that the colder, conductive lid is saltier (and thus denser) than the underlying warm ice. Pits, domes, and small chaos (Michaut & Manga, 2014); ridges (Greenberg et al., 1998); and lobate features (Kattenhorn & Prockter, 2014) have all been suggested to form as the result of cryomagmatism or cryovolcanism. If cryomagmas are sourced from Europa's ocean, which likely has a significant salt content based on the high strength of Europa's induced magnetic field (Kivelson et al., 2000), cryomagmatism could enhance the salt content of the conductive lid. Including a fractional difference in the salt content between the conductive shell and the underlying warm, convecting ice, Δf salt − conv , results in a lower magnitude minimum average density anomaly (Figure 10). Figure 10 shows that the inclusion of Δf salt − conv is well described by a linear combination of the fiducial model and a constant increase in average density anomaly with depth. The magnitude of this increase in average density anomaly with depth is set by the value of Δf salt − conv . This simple behavior holds until ~290 kyr (Figure 10). At this time, the top of the modeled parcel reaches a depth equal to h cond and the difference in salt content of the parcel and the reference material can no longer increase with depth (i.e., Δf salt = Δf salt − conv for every cell in the slab). For ϕ 0 = 0.1, Δf salt − conv > 0.22 is required for the average density anomaly to remain positive at all times meaning the slab would be nonbuoyant (Figure 10). The composition of Europa's ice shell is relatively uncertain, with models of formation and evolution suggesting that the ice shell could be nearly pure water ice or have a significant fraction of hydrated salts (e.g., Kargel, 2000). Although high, a salt content of 22% is conceivable (e.g., Kargel, 2000), and in such a case subduction would be geophysically feasible. Figure 10 Open in figure viewer PowerPoint Effect of salt content, Δf salt − conv , on the time evolution of the average density anomaly in the modeled parcel of subducting material. Assumed fractional differences in the salt content of the conducting and convecting portion of the ice shell, Δf salt − conv , in both the subducting slab and reference shell are indicated by the legend (e.g., Δf salt − conv = 0.1 means that there is 10% more salt in the subducting slab than the underlying convecting ice). All other parameters are unchanged from our fiducial model. So far, we have considered porosity, temperature, and salt content variations with depth and assumed that the ice shell is laterally homogeneous. Assuming that the ice shell is laterally homogeneous leads to initial conditions of a parcel with zero average density anomaly. Although variability in salt content at the optical surface seen in the NIMS data does not directly translate to salt content heterogeneity at depth, this observation acts as motivation to consider laterally heterogeneous salt abundance. Compositional heterogeneity may be produced by thermal and compositional diapirism (Pappalardo et al., 1998; Schmidt et al., 2011), intrusive cryomagmatism (Greenberg et al., 1998; Michaut & Manga, 2014), or cryovolcanism (Kattenhorn & Prockter, 2014). Spatially random emplacement of cryomagma would tend to increase the salt content and density of Europa's ice shell over time, making older regions of the ice shell denser. To explore the effect of lateral compositional heterogeneity, we include a fractional difference between the salt content of the subducting slab and the conductive portion of the average reference ice shell, Δf salt − ref . A Δf salt − ref = 0.01, for example, means that the subducting slab is 1% saltier than the average ice shell and initially ~5 kg/m3 denser than the surrounding material (equation 10). This means that the initial state of our model is one where the parcel is nonbuoyant. When material in the slab is deeper than h cond , the density anomaly from salt is Δf salt = Δf salt − conv , the same as in models shown in Figure 10, which do not include lateral heterogeneity. For material above h cond , however, Δf salt = Δf salt − ref (for the models in Figure 10 this material has Δf salt = 0). Thus, when Δf salt − conv is kept the same, changes to Δf salt − ref only change the average density anomalies out to a time of ~290 kyr. Since all porosity is crushed out by 290 kyr, a comparison of the red curve from Figure 11a and the purple curve from Figure 10 (both corresponding to Δf salt − conv = 0.05) demonstrate that models with different Δf salt − ref have the same average density anomaly after 290 kyr. Figure 11 Open in figure viewer PowerPoint Time evolution of (a) the average density anomaly in the modeled parcel of subducting material and (b) time‐averaged density anomaly of the subducting slab. Combinations of difference in salt content between the conductive slab and warm convecting ice, Δf salt − conv , difference in salt content between the conductive slab and the conductive part of the reference shell, Δf salt − ref , initial porosity, ϕ 0 , and viscosity at the base of the conductive shell, η b , are plotted as indicated by the legend. All other parameters are unchanged from our fiducial model. Δf salt − ref > 0 results in a positive average density anomaly at time zero and is initially nonbuoyant, we define a time‐averaged density anomaly to assess whether subduction is feasible. The time‐averaged density anomaly at a given time, t n = nΔt , is given by (11) Δρ av (t i ) is the spatially averaged density anomaly used throughout this paper at a given time step t i = iΔt and Δt is our time step of 1 year. We define models where subduction is feasible as those with a time‐averaged density anomaly that never becomes negative at any time t n . Because the position and depth of the slab evolve with time (Figures t = 0 ) to the position and depth at time t n . This ensures that from the initiation of subduction, t = 0 , to a time of interest t n , even if some portions of the slab are buoyant, the slab as a whole is always over‐dense and tends to sink. Including lateral heterogeneous salt content allows for slabs that are nonbuoyant. For example, a slab with an initial porosity ϕ 0 = 0.05 , salt content of 2% ( Δf salt − conv = 0.02) , and 1% more salt than the reference conductive lid ( Δf salt − ref = 0.01 ) will be nonbuoyant. Lower η b and ϕ 0 result in higher time‐averaged density anomalies making subduction more likely. Higher Δf salt − conv and Δf salt − ref also result in higher time‐averaged density anomalies. Figure ϕ 0 = 0.1 , η b = 1013 Pa s, Δf salt − conv = 0.05 , and Δf salt − ref = 0.025 the time‐averaged density anomaly remains positive and subduction is possible. Becauseresults in a positive average density anomaly at time zero and is initially nonbuoyant, we define a time‐averaged density anomaly to assess whether subduction is feasible. The time‐averaged density anomaly at a given time,, is given bywhereis the spatially averaged density anomaly used throughout this paper at a given time stepandis our time step of 1 year. We define models where subduction is feasible as those with a time‐averaged density anomaly that never becomes negative at any time. Because the position and depth of the slab evolve with time (Figures 2 and 3 c) this time average is the same as a spatial average over the slab from the location where subduction initiates () to the position and depth at time. This ensures that from the initiation of subduction,, to a time of interest, even if some portions of the slab are buoyant, the slab as a whole is always over‐dense and tends to sink. Including lateral heterogeneous salt content allows for slabs that are nonbuoyant. For example, a slab with an initial porosity, salt content of 2% (, and 1% more salt than the reference conductive lid () will be nonbuoyant. Lowerandresult in higher time‐averaged density anomalies making subduction more likely. Higherandalso result in higher time‐averaged density anomalies. Figure 11 shows a range of parameters for which the time‐averaged density anomaly remains positive at all times and subduction is possible. For a plausible combination of parameters likePa s,, andthe time‐averaged density anomaly remains positive and subduction is possible. In Figure 11 we show models where the modeled parcel of material is buoyant at some times but the slab is always nonbuoyant. It is also possible to construct simulations with plausible parameters in which our modeled parcel of material never becomes buoyant (Figure 12). For example, when ϕ 0 = 0.05, η b = 1013 Pa s, Δf salt − conv = 0.04, and Δf salt − ref = 0.02, the parcel is always nonbuoyant and subduction is feasible regardless of arguments about the buoyancy of the slab as whole (i.e., time averages used above). For larger salt contents and lower porosities than those shown in Figure 12, the modeled parcel will be even more dense and more likely to sink. Figure 12 Open in figure viewer PowerPoint Time evolution of the average density anomaly in the modeled parcel of subducting material showing conditions where the parcel is always nonbuoyant. Combinations of difference in salt content between the conductive slab and warm convecting ice, Δf salt − conv , difference in salt content between the conductive slab and the conductive part of the reference shell, Δf salt − ref , initial porosity, ϕ 0 , and viscosity at the base of the conductive shell, η b , are plotted as indicated by the legend. All other parameters are unchanged from our fiducial model.

4 Discussion In cases where the slab is nonbuoyant, we suggest that the plate will be subducted rather than subsumed (i.e., incorporated into surrounding material as the thermal anomaly diffuses). Even after a thermal anomaly is erased, most models in which subduction is feasible (Figures 11 and 12) have a density anomaly associated with salt content. Unless salts can drain away from the slab on a time scale of ~1 Myr, this material will continue to descend. Moreover, convection models with overturning lids have high surface velocities exceeding 100 mm/yr (Showman & Han, 2005). In these models, lid material descends to the base of the ice shell before the thermal anomaly can diffuse and the scenario is more consistent with subduction than subsumption. Our model greatly oversimplifies the process of subduction. Ideally, we would model the evolution of porosity, temperature, and salt content in a dynamic viscoelastic model or in a convection model that includes plastic yielding (Showman & Han, 2005). In such a model, the evolution of the slab would affect the slab dynamics, including slab bending and evolution of the dip angle. In our model, however, we enforce a constant rate of bending and subduction even if the slab is buoyant. A dynamic viscoelastic model could also include a more sophisticated strain rate and temperature‐dependent rheology for ice (Durham, Kirby, & Stern, 1997; Durham & Stern, 2001; Goldsby & Kohlstedt, 2001). This type of model could also track the evolution of temperature without enforcing strict, but reasonable, temperature boundary conditions at the top and bottom of the slab. These proposed simulations would be computationally expensive, and exploring a large parameter space would be more difficult than with our simple one‐dimensional model. Considering that the evolution of porosity is strongly dependent on temperature, but only weakly dependent on overburden pressure (equations 4 and 5), we do not expect the dynamics of plate bending and subduction to have a strong effect on the buoyancy of the subducting slab. Regardless of these simplifications, our modeling shows that porosity and salt content ultimately determine whether subduction in the ice shell is feasible. We argue that porosity and salt content should be considered in any model of ice shell overturn and subduction on Europa or elsewhere (e.g., O'Neill & Nimmo, 2010; Showman & Han, 2005). Another important issue is that the initial porosity and composition in Europa's ice shell are largely unconstrained. Our modeling shows that relatively small lateral variations in salt content (Δf salt − ref = 0.025) make subduction in Europa's ice shell possible even when porosity is relatively high (ϕ 0 = 0.1) and differences in salt content between the conducting portion of the ice shell and underlying warm ice are reasonable (Δf salt − conv = 0.05). Increasing lateral variation in salt content, overall salt content of the conductive lid, and decreasing porosity all increase the average density anomaly of the slab and make subduction more likely. We note that although we did not explore the effects directly, lateral variations in temperature and porosity would have a similar effect to lateral variations in salt content, producing conditions more favorable to subduction. Possible mechanisms for driving plate motions and initiating subduction are convection (Showman & Han, 2005) and nonsynchronous rotation (Geissler et al., 1998). Our model is independent of these mechanisms and if we assume that subduction is occurring, or has occurred in the recent past, our results provide useful constraints on the composition and structure of Europa's ice shell. Subduction in the ice shell can account for the abundance of extensional features, a lack of compressional features, and young age of Europa's surface (Kattenhorn & Prockter, 2014). Demonstrating that subduction is geophysically feasible under certain conditions allows for further consideration of this exciting possibility: something like plate tectonics on a body other than Earth (subduction likely occurs on Venus, but it lacks lithospheric plates and mid‐ocean spreading ridge analogs (Davaille, Smrekar, & Tomlinson, 2017; Sandwell & Schubert, 1992)). Exchange of surface materials with the subsurface ocean helps determine Europa's astrobiological potential (Chyba & Phillips, 2001; Gaidos, Nealson, & Kirschvink, 1999). Subduction in Europa's ice shell would enhance the rate at which surface material can be delivered to its subsurface ocean, possibly producing conditions that are more conducive to life (Kattenhorn & Prockter, 2014). A better understanding of how salts and porosity are distributed in the ice shell may also shed light on other important processes at work on Europa including diapirism and cryomagmatism (Michaut & Manga, 2014; Pappalardo et al., 1998; Schmidt et al., 2011).

5 Conclusions We have explored the geophysical feasibility of subduction in Europa's ice shell using a one‐dimensional finite difference model that tracks the evolution of temperature and porosity within a subducting slab. We tested the effects of porosity, salt content, conductive lid thickness, basal viscosity, spreading rate, bending radius, and dip angle on the buoyancy of a subducting slab. We found that the buoyancy of a subducting ice slab is most sensitive to salt content and initial porosity. Higher salt contents, lateral variations in salt content, and lower porosity result in slabs that are less buoyant. We find that slabs with salt contents of 5% (Δf salt − conv = 0.05) and 2.5% more salt than the reference conductive lids (Δf salt − ref = 0.025) are nonbuoyant even when they have relatively high initial porosities of ϕ 0 = 0.1. Our work shows that subduction is geophysically feasible under plausible conditions and supports Kattenhorn and Prockter's (2014) geology‐based hypothesis of subduction in Europa's ice shell. Because compositional density variations can remain even after the slab has reached the same temperature as warm ice, material will continue to descend and may transport surface materials to Europa's subsurface ocean.

Acknowledgments We thank S. Kattenhorn and an anonymous reviewer for their helpful reviews. We thank Francis Nimmo for fruitful discussions about this work. The model output and source code needed to reproduce these results are available from the authors upon request (brandon_johnson@brown.edu). The source code needed to replicate our results and model output presented in figures are available through the Harvard Dataverse . https://doi.org/10.7910/DVN/O4A1WU.