On the More Complex Wavelength Dependency of Airy Isostasy in Icy Shells of Ocean Worlds

The topography of ocean worlds is often used to infer ice shell thicknesses by assuming topography is compensated by a basal root. We systematically test the stability of isostatically compensated topography in ice shells. At short horizontal wavelengths, lithospheric strength can support surface topography, while at long wavelengths, buoyancy forces can support topography at the surface and base of the ice shell over geologic time scales. These behaviors are also seen for Airy isostasy in terrestrial worlds. Contrastingly at intermediate scales, the mechanically weak lower ice shell can inhibit the transfer of buoyancy forces to the surface. Factors such as surface temperature can alter the contribution of lithospheric strength, decreasing the stability of a compensating root. This nuanced understanding of icy shell lithospheres provides crucial insights for interpreting surface features and inferring underlying ice shell thickness, with substantial relevance for upcoming space missions to the Jovian system.


Introduction
Topography on terrestrial bodies is generally supported by a combination of factors.These include variations in crustal thickness, known as Airy isostasy (Airy, 1855), variations in density, referred to as Pratt isostasy (Pratt, 1855), and the strength of the lithosphere (Turcotte et al., 1981).The extent to which lithospheric strength contributes to topographic support depends on the wavelength of the topography (e.g., Turcotte & Schubert, 2014).For instance, at short wavelengths, the lithosphere can fully support the topographic load, analogous to a bookshelf easily supporting a single book.However, at sufficiently long wavelengths where the load is spread out across a larger area, the rigidity of the lithosphere is reduced, similar to many books on a sagging bookshelf.In this case, the topography is compensated by isostasy.
The concept of surface topography being compensated by an opposing basal topography (i.e., Airy isostasy) is simple and appealing.It enables interpretation of gravity and topography data collected by spacecraft requiring only an assumption about the density.Consequently, researchers have widely used this concept to estimate variations in the thickness of icy shells on ocean worlds within the Jovian and Saturnian systems.
These icy satellites display a diverse array of topographic features at various horizontal wavelengths.Airy isostasy has been used to estimate ice shell thickness based on observations of surface features and their respective topographies.Shell thickness estimates are commonly derived from long-wavelength topography (e.g., Čadek et al., 2019b;Nimmo et al., 2011;Schenk & McKinnon, 2009), and are often used as an explanation for unexpected topography observations.For example, the unexpectedly large amplitude of Titan's long-wavelength topography may be explained, in part, by isostatic compensation (Nimmo & Bills, 2010).Estimates based on long-wavelength variations in shell thickness assume that the shell is in a constant state of equilibrium.Effectively, ice flow at the base of the ice shell, driven by horizontal pressure differences, reduces basal topography Abstract The topography of ocean worlds is often used to infer ice shell thicknesses by assuming topography is compensated by a basal root.We systematically test the stability of isostatically compensated topography in ice shells.At short horizontal wavelengths, lithospheric strength can support surface topography, while at long wavelengths, buoyancy forces can support topography at the surface and base of the ice shell over geologic time scales.These behaviors are also seen for Airy isostasy in terrestrial worlds.Contrastingly at intermediate scales, the mechanically weak lower ice shell can inhibit the transfer of buoyancy forces to the surface.Factors such as surface temperature can alter the contribution of lithospheric strength, decreasing the stability of a compensating root.This nuanced understanding of icy shell lithospheres provides crucial insights for interpreting surface features and inferring underlying ice shell thickness, with substantial relevance for upcoming space missions to the Jovian system.

Plain Language Summary
Ocean moons have topography on their surfaces that can be used to estimate the thickness of the capping ice shells.In contrast to rocky worlds, the nature of ice shells prevents traditional approaches from being applied.In this study, we use geophysical models to simulate how topography evolves over geologic time.For small features, the stiffness of the ice crust plays a significant role in supporting topography.For larger features, buoyancy forces, like how an iceberg floats, keep topography stable over long periods.In the lower shell near the underlying ocean, the ice is near its melting temperature, which makes it much weaker than the ice nearer the surface.This weakness can prevent the buoyancy forces from reaching the surface for intermediate scale features.Understanding the dynamics of Airy isostasy in ocean worlds can help interpret data collected from space missions, and further our understanding of these icy satellites.and the surface moves in lockstep, maintaining a constant compensation (e.g., Ojakangas & Stevenson, 1989;Stevenson, 2000).
Short-wavelength features are strewn across the icy surfaces of ocean worlds.Europa's surface is most notably characterized by intersecting ridges that can extend for more than 1,000 km with a relief of ∼200 m (e.g., Pappalardo et al., 1999;Prockter & Patterson, 2009).Additionally, dark bands on Europa can rise up to 150 m above the surrounding terrain and span several tens of kilometers (Nimmo et al., 2003b).Although linear features are ubiquitous on Europa, there are also small pits and uplifts with average diameters of around 5 km.Singer et al. (2021) employed a basic Airy isostasy model to estimate a minimum shell thickness for Europa of about 3-8 km, based on the size of these small pits.However, if sufficiently thick, the ice shell likely has the rigidity to support short-wavelength features (e.g., Nimmo et al., 2003a).
Ice shell thickness estimates have also been derived from chaos terrain on Europa.Williams and Greeley (1998) estimated a minimum thickness of just 0.2-3 km at Conamara Chaos, based on the interpretation that the blocks of the chaos terrain were floating in a liquid at the time of formation.In addition to the short-wavelength icy blocks, chaos terrains display longer-wavelength topographic signature suggesting a compensated state (e.g., Schenk & Pappalardo, 2004;Schmidt et al., 2011).
Recently, efforts have been made to refine Airy isostasy for use on ice shells, specifically at low spherical harmonic degrees (see Beuthe, 2021aBeuthe, , 2021b for an extensive review).Hemingway and Matsuyama (2017) argued that the traditional form of Airy isostasy, in which columns of equal width contain equal mass, is not applicable when using spherical geometry as it causes lateral pressure gradients.The authors proposed an approach that assumes continuous pressure along internal equipotential surfaces.Čadek et al. (2019a) assessed the accuracy of these two approaches, along with an "equal stress approach" proposed by Dahlen (1982), to a numerical solution of viscous flow in the crust.The authors concluded that the equal pressure approach may lead to inaccurate estimates of shell thickness for spherical harmonic degrees l ≤ 10 when assuming a constant viscosity.When viscosity of the ice was temperature dependent, all models lost accuracy with increasing spherical harmonic degrees.
A few researchers have presented models that aimed to explain the support mechanism for observed features on the surface of icy ocean worlds such as the plateaus on Titan (Schurmeier et al., 2016) and the megadome on Ganymede (Kay et al., 2018).For both instances, buoyancy forces were not transferred to the surface, which caused the topography at the surface and at the base of the shell to relax away.The megadome on Ganymede was found to be stable over geologic times only when Pratt isostasy was the assumed mechanism of support.On Titan, the plateaus could not be sustained via Airy isostasy, which implies that these features are unlikely to be the result of crustal thickening.
Inherent to the nature of icy shells is the mechanical weakness near the ice-ocean interface.Consequently, a material near its melting point (i.e., the lower ice shell) should be limited in its ability to transfer buoyancy stress to support surface topography.This characteristic injects a level of complexity into the understanding of icy shell topographies and the mechanisms that give rise to their support.In this paper, we explore the mechanics of buoyantly supported topography from short wavelength to hemispheric wavelength scales.With upcoming spacecraft missions such as NASA's Europa Clipper and ESA's JUICE, understanding the dynamics of floating ice shells is crucial for providing a foundation for interpreting the trove of topographic data that could profoundly alter our knowledge of icy ocean worlds.

Finite Element Simulations
We use the commercially available Hexagon Marc finite element package to test systematically for buoyant support of topography in icy shells.The software is well vetted for investigations of lithospheres of icy shells (e.g., Damptz & Dombard, 2011;Dombard & McKinnon, 2000, 2006a;Kay & Dombard, 2018).The code simulates Maxwell viscoelasticity, capturing the behavior of geologic materials at both short and long timescales, (and while not invoked in our simulations, it can include plasticity as a continuum approximation for brittle failure at large stresses).We test for buoyant support of topography for wavelengths λ = 10 km up to hemispherical scales (i.e., spherical harmonic degree l = 2), increasing logarithmically.Although our study is designed to be applicable to all icy ocean worlds, we use Europa to guide the selection of our model parameters.The full finite element simulation is comprised of two steps.A thermal simulation is run to determine the temperature structure, and the results of the thermal simulation are used in the mechanical simulation.

Mesh and Boundary Conditions
We utilize planar half space meshes split into two layers: an outer ice shell, modeled to be 20 km thick based on estimates for Europa (e.g., Billings & Kattenhorn, 2005;Hussmann et al., 2002), and an ocean with a boundary far enough away to not affect the solution (see Figure S1 in Supporting Information S1 for model schematic).For short wavelengths, the ocean mesh layer is ∼2.5 times the wavelength.For long wavelengths (l ≥ 6), the thickness of the ocean layer is reduced to <2.5 times the wavelength, because of computational limits for large meshes.The planar meshes we implement do not include support from membrane stresses.Membrane stresses become appreciable when the horizontal scale of a spherical harmonic degree, l, is roughly the same as the value given by  ≈ √ ∕ , where R is the planetary radius and d is the lithospheric thickness (Turcotte et al., 1981).For Europa, this corresponds to a spherical harmonic degree of ∼16 (∼600 km).As will be discussed in subsequent sections, the additional support provided by membrane stress simply adds to the already sufficient support to maintain topography at those wavelengths.
Elements in our mesh are uniformly distributed in the ice shell layer.In the ocean layer, elements are more densely packed near the ice-ocean boundary, ensuring a higher resolution in this region, while still adhering to a 4:1 aspect ratio throughout.Mesh resolution is variable depending on the scale of the wavelength.At the smallest scales, the ice portion of the mesh has, at minimum, 400 finite elements (40 in depth and 10 in width).At the largest scales, this ratio is reversed to keep within the 4:1 aspect ratio, with 400 elements in width and 10 in depth.Simulations with increased element density in the ice shell portion of the mesh showed no noticeable differences in the results, confirming that the mesh resolution used is adequate for capturing the key dynamics across all scales.Free-slip conditions are applied to the symmetry boundaries at the sides of the mesh, restricting normal displacements, and the bottom of the mesh is locked horizontally and vertically.On the surface, we place sinusoidal topography with an amplitude h of 100 m for all wavelengths, as wavelength controls buoyant support rather than topographic amplitude (see Turcotte & Schubert, 2014).We simulate two scenarios for initial topography at the base of the ice shell: (a) a flat basal surface (i.e., no initial compensating topography), and (b) a root in isostatic equilibrium.The maximum amplitude of the isostatic root H is determined by the density contrast between the water ice and ocean where (1) In this equation, the density of the liquid water ρ w is assumed to be 1,000 kg/m 3 , while the density of the ice shell ρ i is 920 kg/m 3 .Gravitational body forces are applied uniformly to the mesh with an acceleration equal to 1.315 m/s 2 .

Material and Thermal Properties
The elastic parameters for ice include a Young's modulus of 9.32 GPa and a Poisson ratio of 0.325 (Gammon et al., 1983).However, in simulations, applying a gravitational load to a compressible material leads to deviatoric stresses from self-compaction.Following Dombard and McKinnon (2006b), we make the material nearly incompressible by adjusting the Poisson ratio to approach 0.5, while also reducing the Young's modulus to 7.83 GPa as to maintain the flexural rigidity of the ice (Turcotte & Schubert, 2014, Equation 3-72).The viscosity structure of the ice layer has the greatest influence on relaxation rates (Dombard & McKinnon, 2006a), and when incompressibility is approximated, the variation is negligible (Dombard et al., 2007).
Additionally, we follow the ductile creep laws from Goldsby and Kohlstedt (2001).Several mechanisms affect the viscous flow of ice below the brittle-ductile transition.There are three different dislocation (dis) creep mechanisms with temperature dependent contributions (Durham & Stern, 2001), diffusion (diff) creep; and dislocation in an easy-slip (ES) system and grain boundary sliding (GBS) acting as rate limiters.Thus, the total strain rate is Each of these strain rates can take the general form where    is the equivalent strain rate, A is a material dependent pre-exponential constant normalized for uniaxial deformation, m is a grain size index relating to the grain size d, Q is the activation energy, R is the universal gas constant, and T is absolute temperature.We assume a grain size of 1 mm (e.g., Dombard & McKinnon, 2006a).The full creep parameters are listed in Durham et al. (2010) and references therein.
The thermal conductivity of water ice is temperature dependent, given by 651/T (Petrenko & Whitworth, 1999).In our steady state thermal simulation, we lock the surface temperature at 100 K and the ice-ocean interface temperature at the melting point of ice, 273 K (q ≈ 32 mW/m 3 ).Heat flows on Europa might be higher, instigating convection near the base of the shell (e.g., McKinnon, 1999;Tobie et al., 2003).Nonetheless, for our purposes, the situation is mechanically identical: a structurally weak base overlaid by a relatively thinner lithosphere (see below).
The water below the ice shell lacks material strength and therefore does not provide resistance; rather, it provides a buoyant restoring force.The viscosity is so low compared to the lithospheric ice in the upper shell that it simply needs to accommodate space for the deforming ice at least as fast as that ice is deforming.Therefore, we set the material properties of the water to be the same as the overlying ice layer, with the exception of the density contrast.We also test for the effects of surface temperature, which have been shown to be a controlling factor in creep rates (e.g., Damptz & Dombard, 2011), and the effects for a reduction in the melting temperature of ice due to the presence of an antifreeze such as ammonia (e.g., Spohn & Schubert, 2003).To test for effects of buoyancy, we run the simulations without the density contrast, precluding isostasy.Variations in model parameters are shown in Table S1 in Supporting Information S1.
We run the simulations for 3 Gyr, a time more than sufficient to capture any appreciable flow in the lower ice shell.Time steps are controlled by the minimum Maxwell time, which is dependent upon, in part, the viscosity (Turcotte & Schubert, 2014).We limit minimum viscosity in the mesh to 10 18 Pa s due to computational constraints.While a reduced viscosity cutoff accelerates relaxation (Figure S2 in Supporting Information S1), the core behavior remains consistent.
Simulations are performed under full large-strain formulation, which includes a second-order term to strain-displacement relationships (e.g., Ranalli, 1995), and a continual geometric update.Although strains are typically small, the large strain scheme is essential for the geometric update, as without it, stresses would not decay as topography relaxes (e.g., Dombard & McKinnon, 2006a).In addition, we apply constant dilatation across the elements to eliminate potential numerical errors caused by simulating a nearly incompressible material.

Results
The evolution of topography can be described best by three general wavelength categories: short, intermediate, and long.Figure 1 shows the simulation results for the evolution of topography at the surface and base of the ice shell for three representative wavelengths within these categories.
At short wavelengths, the surface topography remains stable over the simulation time regardless of buoyancy.In contrast, the basal topography relaxes away, losing 991 m of its initial amplitude.For intermediate wavelengths, the surface topography relaxes in both the buoyancy and no buoyancy simulations.However, the degree of relaxation is substantially greater when the density contrast is not present.For instance, for a wavelength of 100 km, the surface loses approximately 65% of its elevation over the simulation time, compared to an 85% loss for the no-buoyancy scenario (Figure 1c).The root topography at intermediate wavelengths also relaxes when the density contrast is present.Moreover, the root relaxes at a different rate than the surface topography, leading to a deviation from the expected isostatic equilibrium (Figure 1d).This discrepancy arises from the shared contribution of both lithospheric support and buoyancy (Figure 2a).In the case of long wavelengths, without the density contrast, the surface completely relaxes away in less than 100 Kyr.Conversely with the density contrast, the surface and the root are stable over geologic times.
In Figure 1d, root amplitudes serve as indicators of the varying contribution of lithospheric strength and buoyancy across different wavelengths.At the shortest wavelengths, lithospheric strength dominates, and the root is decoupled from the surface and can flow away.Thus, the predicted root topography, found using Equation 1 on the simulated surface topography, remains close to the initial amplitude of 1,150 m, while the simulated root topography collapses.As wavelengths reach the intermediate scale, lithospheric strength continues to be the primary support mechanism, but buoyancy begins to contribute.This effect is evident around 100 km wavelength, where the simulated and expected root amplitudes are similar, but their rates of amplitude loss differ, indicating that isostatic equilibrium is not maintained.Beyond these scales, where lithospheric strength becomes less effective, the rate and magnitude of amplitude loss for both the simulated and expected root decreases.Ultimately, at the longest wavelengths, buoyancy becomes the sole support mechanism, and the simulated and expected root are near the initial 1,150 m amplitude.
Simulations with a surface temperature of 100 K show that the lithosphere supports the topographic load up to wavelengths of ∼30 km (Figure 2).Lowering the surface temperature increases the wavelengths at which lithospheric strength can effectively contribute to topographic support.A decreased surface temperature also reduces the degree of relaxation occurring at intermediate wavelengths.For example, at 100 km, the simulated surface topography with temperatures of 100 K relax ∼65 m, compared to ∼15 m when the surface temperature is 70 K.A similar effect is observed when the melting temperature (i.e., the temperature at the ice-ocean interface) is reduced to 176 K (Figure 2b).Under these conditions, the contribution of lithospheric support is greater at longer wavelengths, and the maximum topographic relaxation is reduced.All curves originate at time = 0.For each wavelength, the solid line is the root amplitude of the buoyancy simulation, and the dashed line is the expected root amplitude based on the surface topography of the buoyancy simulation.When the dashed line and solid line are the same, the isostatic equilibrium is maintained.

10.1029/2023GL106046
6 of 9 At wavelengths <80 km, the simulations with an initially flat ice shell base are nearly identical to the simulations without the density contrast (Figure 2b).This finding further underscores that lithospheric strength is the primary mechanism of topographic support at short wavelengths.However, at wavelengths greater than 100 km, without an initial compensated root, the surface topography relaxes to just a few meters over the simulated time.At the longest wavelengths, a root forms that is in isostatic equilibrium with the reduced surface topography.

Discussion
The prevalent assumption that Airy isostasy can be utilized for estimating ice shell thickness does not hold true except at the longest wavelengths.A significant oversight in this assumption lies in the positioning of the lithosphere relative to the ice-ocean boundary.On Earth, the lithospheric boundary is situated beneath the density The "No Strength" curve is the difference between "Buoyancy and Strength" and "No Buoyancy."Because the simulations (except the "Flat root" ones) are set up to be in isostatic equilibrium, the ratio should be pegged to a value of 1 if the precepts of Airy isostasy hold; that is, Airy isostasy fails when the amplitude ratio is <1.Refer to Section 2 and Table S1 in Supporting Information S1 for details on the boundary conditions of each simulation type.contrast of the crust-mantle boundary.Thus, when a surface load deforms the lithosphere, both the crust and the upper mantle will respond to the deformation (e.g., Barrell, 1914;Daly, 1940).However, in the case of ice shells, the lithosphere is embedded within the ice itself, because the ice in the lower portion of the shell is approaching its melting point and is therefore weak over geologic time scales.A convective portion in an ice shell that passes a higher heat flux would only exacerbate this effect.Relative to a conductive ice shell, a convective shell of the same total thickness would see higher temperatures at shallower depths, resulting in a thinner lithosphere and a thicker channel that decouples the root topography from the surface, permitting lower crustal flow to longer wavelengths.The thinner lithosphere would shift the "No Buoyancy" (i.e., all lithospheric strength) curves in Figure 2 to the left, while the thicker channel would shift the "No Strength" (i.e., all buoyancy) curves to the right.The combined effects would yield an exacerbated (i.e., wider) intermediate wavelength zone of weakness.
Our findings also illustrate that the lithosphere of ice shells can support topographic loads at short wavelengths.The short-wavelength scenario shows surface topography to be stable with and without buoyancy, while the initially isostatic basal topography relaxes away (Figure 2a), indicating that contribution of buoyancy supporting the load is insignificant and that any root sticking into the ocean would flow away rapidly.The length at which lithospheric support begins to decrease is contingent on the surface temperature, and potentially, the presence of antifreeze in the ocean.These findings agree with Damptz and Dombard (2011), who highlighted the significant role of surface temperature in controlling creep at the base of the lithosphere for the moons of Jupiter and Saturn.Due to viscous creep, the lithosphere is not static.Stresses relax at the base of the lithosphere, allowing the lithosphere to thin with time.Surface temperature, and to a lesser degree heat flow and grain size, control creep rates.
It is important to note that while the surface topography at both short and long wavelengths can be stable over geologic times, the basal topography exhibits distinctly different behaviors at varying spatial scales.At short wavelengths, basal topography relaxes away indicating topographic features at these scales are not compensated.Conversely, at long wavelengths, topography at the base of the ice shell stays in isostatic equilibrium with the surface topography.This finding can be attributed to the relative size of the vertical flow region and the lateral extent of the topography.For long wavelengths, the vertical extent of the area where material can flow laterally and redistribute is small compared to the horizontal extent (cf.Stevenson, 2000).This imbalance results in a highly channelized flow that inhibits the overall relaxation process observed at shorter wavelengths.
Moreover, when the wavelength of a load is small compared to the planetary radius, the stresses caused by the body's curvature can be ignored.However, when curvature is introduced, horizontal membrane stresses can contribute to the support of the surface load (Turcotte et al., 1981).In the context of this study, membrane stresses are appreciable at approximately spherical harmonic degree l = 16, assuming the lithosphere is ∼6 km thick.This degree corresponds to a wavelength of ∼600 km.At this scale, topography loses just 10% of its initial height when surface temperature is 100 K, and this relaxation is further reduced for colder surface temperatures (Figure 2).Consequently, the inclusion of membrane stresses would reinforce our findings regarding the stability of long wavelength topography.
To sum up, there are two competing effects for the application of Airy isostasy in an ice shell over an ocean, one of which does not exist in traditional applications of Airy isostasy in terrestrial worlds.First, the lithospheric strength of icy shells is a sufficient support mechanism for topography at short wavelengths, without needing help from buoyancy.This finding is evidenced by the "No Buoyancy" curve in Figure 2a, and it mirrors the lithospheric behavior of rocky planets.In contrast, the "No Strength" curve in Figure 2a, which does not come from simulations but from subtracting the "Buoyancy and Strength" and "No Buoyancy" curves and thus mimics a scenario where any lithosphere would have no strength, stands in contrast to the expected behavior in a terrestrial world, which would be locked at 1 at all scales.These two curves-"No Buoyancy" and "No Strength"-together determine the final topographic support.
To understand this behavior, consider the effect of pressing one's finger into a seat cushion.The cushion diffuses the stresses introduced by the force of the finger, and the opposite side of the cushion is negligibly deformed, if at all.A sufficiently wide scale load, however, would feel the hard seat underneath, even when pressed into the cushion at the same force per unit area.Translating to ocean worlds, the seat cushion inhibits buoyant support at short wavelengths, though the lithosphere can provide the support, while concurrently, the basal topography can relax away without affecting the surface.At long wavelengths, the seat cushion is too thin for buoyancy to be inhibited, so topography can be supported, even if the lithosphere provides no support.Simultaneously at these scales, highly channelized flow precludes relaxation of the basal topography.At intermediate scales, neither the seat cushion or the lithosphere are fully effective, and topography cannot be supported over geologic timescales by a crustal root alone, and would need another mechanism of support (e.g., Pratt).How these curves add, and thus the size of the intermediate zone, is dependent on local conditions, as suggested by Figure 2b.For the case of Europa explored here, this wavelength range spans over an order of magnitude (from <100 km to >1,000 km), with substantial lack of support for equatorial regions where surface temperatures are higher.
This study begins to address the complexities of topographic support mechanisms in icy shells.The results, particularly those regarding varied surface and ocean temperatures and rheological conditions, offer insights into how other icy bodies might respond, and serve as a foundational point for more extensive future research that could delve into a broader range of parameters, including different assumptions about viscosity and rheology.Additionally, thought experiments on how different conditions would affect the lithospheric thickness and the thickness of the decoupling flow channel, and thus the impact on the curves in Figure 2 (such as was done for the above discussion of a convective ice shell), could also provide insight on how the phenomenon would play out on different ocean worlds.

Conclusion
Our results yield a more nuanced understanding of the wavelength dependence of icy shell evolution, offering insights that could enhance the interpretation of topography and gravity data from the icy ocean worlds.While surface topography of short-wavelength features persists over geologic time, compensated roots relax away quickly.Consequently, shell thickness estimates derived from short-wavelength features will not accurately represent the local thickness.In contrast, long-wavelength topography could be supported by shell thickness variations, such as on Titan (Nimmo & Bills, 2010).Surface topography at intermediate scales is unlikely to be isostatically compensated via buoyant support.
The practical implications of this more complex behavior have already been revealed at other icy worlds.The isolated plateaus on Titan appear to exist in the short wavelength regime, being supported by the strength of Titan's lithosphere with any root relaxing away quickly, suggesting top-down construction of the plateaus (Schurmeier et al., 2016).Conversely, the megadome on Ganymede appears to be firmly in the intermediate regime, unable to be supported by either lithospheric strength of a buoyant crustal root and implicating Pratt isostasy (Kay et al., 2018).
Such complexities and limitations underscore the challenges inherent in applying traditional isostatic models to icy shells and highlight the importance of continuous model refinement.Furthermore, the importance of a nuanced understanding of the dynamics of icy shells cannot be overstated, especially in the context of future space missions.

Figure 1 .
Figure 1.(a) 30 km, (b) 100 km, and (c) 1,635 km (spherical harmonic degree l ≈ 6) simulations of surface topography.The initial topography (black dashed line) is consistent across each model with sinusoidal amplitudes of 100 m.The left side of x = 0 are simulations with an ice-ocean density contrast (i.e., "Buoyancy"), while the right side are simulations without the density contrast (i.e., "No Buoyancy").The progression of time is indicated by the shade of the color, with lighter shades representing earlier time steps and darker shades indicating later ones.(d) Time series of root amplitudes for the wavelengths presented in (a)-(c).The color scheme for the amplitudes matches that of the wavelengths in panels (a)-(c).Additionally, root amplitudes for 55 and 316 km are shown in dark gray and light gray, respectively.All curves originate at time = 0.For each wavelength, the solid line is the root amplitude of the buoyancy simulation, and the dashed line is the expected root amplitude based on the surface topography of the buoyancy simulation.When the dashed line and solid line are the same, the isostatic equilibrium is maintained.

Figure 2 .
Figure2.Ratio of initial to final topographic amplitude for simulations across simulated wavelengths for (a) simulations with and without buoyancy and (b) variations in the assumed initial state of the temperature at the surface and base of the ice shell.The "No Strength" curve is the difference between "Buoyancy and Strength" and "No Buoyancy."Because the simulations (except the "Flat root" ones) are set up to be in isostatic equilibrium, the ratio should be pegged to a value of 1 if the precepts of Airy isostasy hold; that is, Airy isostasy fails when the amplitude ratio is <1.Refer to Section 2 and TableS1in Supporting Information S1 for details on the boundary conditions of each simulation type.