Elastic-plastic flexure and heat flux on icy moons

Abstract

Investigating ice shell flexure observed in elevation data provides one of the few ways to infer the ice shell structure and the thermal evolution of icy moons. Most previous work have relied on the classical elastic plate model to fit observations and infer the elastic thickness of the lithosphere and heat flux at the time of deformation. Here, we show how to compute more realistic elastic-plastic flexural profiles considering ice friction in the upper ice shell and plastic flow at deeper levels. This approach does not rely on the elastic thickness mathematical construct and can be used to directly infer the ice shell heat flux and/or rheology from observed flexural profiles. We show that, under high plate curvature (>10⁻⁵ m⁻¹), purely elastic models predict unrealistic oscillations near and in the flexural bulge region, which can bias the estimated elastic thickness and heat flux. Applying our elastic-plastic flexural model to Ariel and Ganymede, we find heat fluxes lower than previous estimates, which has implications for our understanding of the geodynamic history of these worlds. We demonstrate that the maximum curvature of the synthetic elastic plate profile should be used when relating elastic thickness to heat flux, and discuss heat flux estimation on icy moons.

1 Introduction

Investigations of ice shell flexure, as observed from stereo-derived topographic profiles, have been commonly used to provide information on the interior structure and evolution of icy moons, including Ganymede (e.g., Nimmo & Pappalardo, 2004; Nimmo et al., 2002), Europa (e.g., Billings & Kattenhorn, 2005; Hurford et al., 2005), Enceladus (e.g., Giese et al., 2008), Dione (e.g., Hammond et al., 2013), and Ariel (e.g., Beddingfield et al., 2022; Peterson et al., 2015). The most commonly used approach is to fit the observed flexure to an elastic plate model to infer the effective elastic thickness (hereafter denoted as elastic thickness), which provides information on the rigidity of the body’s ice shell. The widespread use of this approach lies in its simple analytical expression (e.g., Turcotte & Schubert, 2002), allowing to explore a wide parameter range at multiple locations. With improved topography models for icy satellites (e.g., Enya et al., 2022; Park et al., 2024; Schenk & McKinnon, 2024), flexural analyses can be applied across the surface to infer lateral variations in the interior structure of the body. More sophisticated models, accounting for time-dependent elastic–viscous–plastic rheology have also been proposed (e.g., Damptz & Dombard, 2011). However, while these models have shown the importance of the time-dependent response of the lithosphere, with implications for the interpretation of the inferred elastic thickness, they require using non-trivial finite-element solvers (see also Dombard & McKinnon, 2006).

For geologic interpretations, the elastic thickness can be converted to a heat flux using several approaches. For example, this can be achieved by setting the bending moment of the elastic plate equal to the bending moment of a more realistic plate with a temperature-dependent rheology that considers fracturing and plastic flow (e.g., Burov & Diament, 1995; McNutt, 1984). One critical issue with this approach is related to the seemingly arbitrary selection of the input plate curvature, which affects the calculation of the bending moment (Goetze & Evans, 1979). The curvature of the plate varies with position by orders of magnitude and has a prominent effect on the inferred bending moment and associated heat flux. Some previous works have used the curvature at the first zero-crossing of the flexural profile (e.g., Smrekar et al., 2022) or the maximum plate curvature (e.g., Nimmo & Pappalardo, 2004), but the reliability of such assumptions has not been investigated in the framework of ice shell deformation (see also Mueller & Phillips, 1995). Alternative approaches have assumed the base of the elastic lithosphere to be defined by a rheology-dependent isotherm in combination with a specific Deborah number (Hammond et al., 2013; Nimmo et al., 2002). The Deborah number is defined as the ratio of the stress magnitude (viscosity times strain rate) and the elastic shear modulus, and this parameter weighs the elastic contribution to the plate’s rheology as a function of the loading time scale (e.g., Mancktelow, 1999). However, it remains unclear what Deborah number should be assumed when the plate is elasto-plastic (rather than purely elastic) and whether both the equal-bending moment and the Deborah number approaches lead to similar results.

In addition, for most icy moons, the retrieved elastic thicknesses were generally found to be small (\(<\)5 km) indicating high heat fluxes (\(>\)100 mW m\(^{-2}\)) at the time of the deformation (e.g., Nimmo, 2018; Van Hoolst et al., 2024). Such information are critical to assess icy moons’ habitability (e.g., Grasset et al., 2013; Vance et al., 2023) or the dynamics of interior plumes (Berne et al., 2024; Spitale et al., 2025). However, one key question is whether the elastic plate model can be used to reliably predict the flexure of an ice shell (and the associated heat flux) with a thin elastic lithosphere, wherein plasticity might dominate. Finally, the presence of impurities or water ice layers (e.g., Buffo et al., 2020; Souček et al., 2023) will affect the yield strength and flexural response of the ice shell, but such features cannot be included onto simple analytical elastic plate models.

In this work, we follow the framework developed in Mueller and Phillips (1995), to test the applicability of elastic plate models on icy satellites. We compare elastic plate models to elastic-plastic flexure solutions using an example yield strength envelope for Europa and provide a criterion for curvature selection when relating elastic thickness to heat flux. We discuss synthetic plastic-elastic flexure profiles considering ice shell impurities and demonstrate that local weak zones (such as from liquid pockets) affect the ice shell’s flexural response. We further show that previous work that used a Deborah number overestimated the heat flux, with examples for Ariel (Peterson et al., 2015) and Ganymede (Nimmo et al., 2002).

2 Methods: Plate flexure

The deflection of a plate can be described using the von Kármán equation (Fung, 1965). Under the assumption of negligible shear force applied to the top and bottom of the plate, which is adequate for icy satellites, the von Kármán equation can be simplified as (e.g., Mueller & Phillips, 1995; Turcotte & Schubert, 2002)

\begin{equation} \frac {d^2M(x)}{dx^2} = -N(x) \frac {d^2w(x)}{dx^2} + \Delta \rho gw(x) + \sigma _{zz}(x,0). \tag {1}\label {eq:Eq1} \end{equation}

In this equation, \(x\) is the horizontal dimension, \(M\) is the bending moment of the plate, \(N\) is the in-plane force, \(w\) is the deflection, and \(\sigma _{zz}\) is the vertical normal stress applied to the plate. Other parameters include the density contrast of materials below and above the plate, \(\Delta \rho \), and the local gravitational attraction, \(g\). We note that the model is invariant with respect to the other horizontal direction. More sophisticated 2D elasto-plastic flexure modeling is out of the scope of this study.

Under elastic flexure theory, the bending moment varies linearly with plate curvature and is given by the product of flexural rigidity and curvature. A boundary condition is also typically applied such that flexure at \(x=0\) equals to some constant \(w_0\) that is determined from observations. For a broken plate, the curvature (or bending moment) also equals zero at \(x=0\). Ignoring in-plane forces, this leads the above equation to only depend on a flexural parameter \(\alpha \) and the central deflection \(w_0\) as (see Turcotte & Schubert, 2002).

\begin{equation} w(x)=w_0 e^{-x/\alpha } \cos {} (x/\alpha ) \tag {2}\label {eq:Eq2} \end{equation}

where the flexural parameter can be expressed as

\begin{equation} \alpha = \left (\frac {ET_e^3}{3(1-\nu )\Delta \rho g}\right )^{1/4}. \tag {3}\label {eq:Eq3} \end{equation}

In this equation, \(T_e\) is the elastic thickness of the lithosphere, \(E\) is Young’s modulus, and \(\nu \) is Poisson’s ratio. The above assumptions have been widely used to infer the strength and thickness of the elastic lithosphere when and where the gravity field were insufficiently known, including on Earth (e.g., McNutt & Menard, 1982), Venus (e.g., Smrekar et al., 2022), Mars (e.g., Grott et al., 2005), and icy moons (e.g., Nimmo et al., 2002).

For elastic-plastic flexure, however, the moment–curvature relationship is non-linear and solutions to Equation (1) need to be determined numerically. Following Mueller and Phillips (1995), we first determine the moment–curvature relationship based on an input interior yield strength envelope (YSE). The bending moment is integrated considering both tensile and compressional yield stresses together with a curvature and associated neutral plane depth (e.g., Goetze & Evans, 1979). Bending moments are estimated for a wide range of curvatures, and the non-linear relationship between the two is parametrized using spline functions (Woltring, 1986). We note that throughout the rest of the manuscript we ensure that all constructed moment–curvature spline functions adequately reflect the input data. Finally, Equation (1) is numerically solved as a two-value boundary problem using a substitution of variables approach (see Phillips, 1990). Codes to solve the elastic-plastic flexure equation for an input YSE, but also to convert an input YSE to \(T_e\) (and reciprocally), are archived in Broquet (2026) and Broquet (2022).

3 Results

3.1 Curvature selection when relating yield strength envelopeand elastic thickness

In this section, we compare flexural profiles of a broken elastic plate (Equation (2)) to our YSE-dependent elastic-plastic flexure model. Here, we use Europa’s ice shell as an example considering a single central flexure (\(w_0\)) and assumed heat flux, but note that a wider parameter search is discussed further below. For Europa, we use the YSE described in Dombard and McKinnon (2006), which considers Byerlee’s law for ice in the upper brittle ice shell as well as dislocation and grain boundary sliding in the lower ice shell (Figure 1a). Assumed parameters include a surface temperature (\(T_s\)) of 80 K, a thermal gradient of 18 K/km, a strain rate of 10\(^{-15}\) s\(^{-1}\), and a 1 mm grain-size. All relevant rheological and interior parameters are summarized in Table S1. For this YSE, the associated moment–curvature relationship is seen to saturate for curvatures greater than |3.5|\(\times \)10\(^{-6 }\) m\(^{-1}\) (Figure 1b). After this threshold limit, yielding occurs and any additional increase in the plate curvature does not induce a significant increase in the bending moment that is supported within the plate. The saturation threshold depends strongly on the yield strength envelope, but can be estimated by simply estimating the total bending moment of the plate.

Figure 1: Example YSE for Europa (a) together with its associated moment–curvature relationship (b). For comparison, we show the analytical elastic moment–curvature relationship considering a 3 km elastic thickness. The YSE is from Dombard and McKinnon (2006).

As expected, this saturation behavior is not reproduced by the analytical elastic curvature-moment relationship, which allows the elastic plate to maintain unphysically high curvatures (Figure 1b). Interestingly, we further observe the elastic and elasto-plastic models to slightly differ in terms of bending moment before the saturation point. Although plate flexure theory shows that elasto-plastic plates have an elastic core located near the midpoint of the plate (e.g., Burov & Diament, 1995), even for small amounts of yielding, a plate will not behave fully elastically. Plasticity near the surface and at the base of the plate will both affect the bending moment. We note that Europa’s YSE and associated plate saturation point strongly depend upon the assumed ice shell rheology and temperature. A full exploration of such parameters is out of the scope of this study.

The constructed YSE and associated curvature-moment relationships are used to compute elastic-plastic flexure (see Equation (1)). Assumed parameters include a central deflection, \(w_0\), of 0.5 km, Poisson’s ratio of 0.33, Young’s modulus of 1 GPa, a mean gravitational attraction of 1.314 m s\(^{-2}\), and a density contrast of 1000 kg m\(^{-3}\) (Table S1). Because we are comparing synthetic data, these assumptions do not affect the main results presented here. The estimated elastic-plastic flexure solution is then compared to a broken elastic plate model, for which one needs to assume an elastic thickness. To estimate the elastic thickness, we use the classical YSE to elastic thickness formalism (e.g., McNutt, 1984). Therein, the curvature-dependent bending moment of the YSE is equated to an elastic bending moment. As the elastic bending moment analytically depends on the flexural rigidity and plate curvature (assumed to be constant), this approach can be used to link a YSE and assumed plate curvature to an elastic thickness, or vice versa. However, and as discussed above, a key unknown is what curvature one should use to adequately infer the elastic thickness. Here we consider a wide range of curvatures to estimate several elastic thicknesses using the above formalism. Elastic thicknesses are then used to construct synthetic broken elastic plate models using the same parameters as for the elastic-plastic solutions. A root-mean-square (rms) misfit is then used to compare broken elastic plate and elastic-plastic flexural profiles. The curvature used to construct the best-fit broken elastic plate model is retained and is then compared to the laterally variable curvature of the synthetic profiles. This approach enables the inference of the representative curvature required to relate the elastic thickness of the lithosphere to a yield strength envelope (and heat flux), and conversely, to derive an elastic thickness from a yield strength envelope and input flexural parameters.

Figure 2: Elastic-plastic versus broken elastic plate flexure for Europa. (a) Root-mean-square misfit between the elastic-plastic and synthetic broken elastic plate models for various elastic thickness and curvatures. (b) Flexural models for broken elastic (orange) and elastic-plastic (black) plates. The orange shade shows all tested broken elastic plate models. (c) Absolute value of the curvature as a function of distance for the best-fit models. The curvature leading to the lowest misfit between the broken elastic and elastic-plastic plate is shown by the red horizontal line. The x-axis in (b, c) corresponds to the distance away from the load.

For curvatures ranging from 10\(^{-7}\) to 10\(^{-5}\) m\(^{-1}\), the elastic thickness associated with our example YSE ranges from 5.0 to 2.3 km (Figure 2a). The broken elastic plate model that best-fits the elastic-plastic model has an elastic thickness of 3.0 km, and was retrieved for a curvature of 3.6\(\times \)10\(^{-6}\) m\(^{-1}\) (Figure 2b). This optimal curvature is best estimated considering the maximum curvature of the best-fit broken elastic plate model (Figure 2b,c), in agreement with previous work (e.g., Mueller & Phillips, 1995). As expected, such curvature is found near the midpoint between the flexural profile start and the first zero-crossing flexure point. If instead, the curvature were selected at the first zero-crossing point (as in Smrekar et al., 2022), here with a value of 2.3\(\times \)10\(^{-6}\) m\(^{-1}\), the elastic thickness of the lithosphere would be overestimated by 10%.

We have also explored a wider parameter space considering positive and negative \(w_0\) ranging from -4 km to 4 km as well as heat fluxes of 20 to 300 mW m\(^{-2}\), building on the Ariel and Ganymede investigations presented below (Figure S1). In all tested scenarios, we look for the input curvature that allows for a best-fit between the elastic model (which depends on an elastic thickness and curvature) and the elastic-plastic solution (which depends on a heat flux). The best-fit curvature is then compared to the maximum curvature of the synthetic elastic model, as was done above. This analysis shows that for a wide range of heat fluxes (thus yield strengths), as well as central deflections (thus curvatures), selecting the maximum curvature of the synthetic elastic profile is always a good approach to match elastic-plastic and elastic models and therefore to estimate the heat flux from an elastic thickness.

These analyses demonstrate that a temperature-dependent yield strength envelope can be related to an elastic thickness by considering the maximum curvature of the synthetic profile that best fits an observed flexural profile. In the tested scenario for Europa with a relatively small central deflection of 0.5 km, flexure obtained by the broken elastic plate model is also found to closely resemble that from the elastic-plastic model.

3.2 Non-elastic flexure

As discussed above and shown in Figure 1b, the bending moment of a lithosphere that considers fracturing and plastic flow will start to saturate under high curvatures. Where yielding occurs, an elastic-plastic flexural profile will differ from the elastic plate model that can sustain infinitively high curvatures. To illustrate this effect, we have increased \(w_0\) to 1 and 2 km, as greater central deflection leads to increasing plate curvatures (Figure 3a,b). We note that such high central deflection have been observed on Ariel, with values reaching \(\sim \)4 km (Peterson et al., 2015). As in the previous section, the maximum curvature in the synthetic profile should again be used to determine the best-fit elastic thickness associated with the YSE and input deformation parameters (Figure 3c,d).

Figure 3: Effect of plasticity in the ice shell. Elastic-plastic versus broken elastic plate flexure for Europa under different maximum central flexure \(w_0\) of 1 km (a) and 2 km (b). The associated curvatures are shown in (c) and (d) with the best-fit shown as the red dashed line and the YSE saturation point displayed as the solid grey line. The figure uses the same format and color-scheme as Figure 2. The x-axis in all panels corresponds to the distance away from the load on either side.

We observe that yielding of the plate before the first zero-crossing point and in the flexural bulge region both lead to a flexural profile that strongly departs from the broken elastic plate model (Figure 3a,b). A large flexural bulge with positive relief and high curvature, above the theoretical saturation point (Figure 3c,d), is seen in the elastic model, but absent from the more realistic elastic-plastic solution. This effect is more pronounced as \(w_0\) increases.

This analysis shows that under high flexural curvatures, elastic and elastic-plastic solutions can substantially differ. In particular, by neglecting plasiticity, elastic flexure modeling will lead to the apparition of unrealistic plate oscillations in the flexural bulge region. Thus, the fit of observed data in such region, which has been typically investigated by prior elastic flexure models (e.g., Hurford et al., 2005; Nimmo et al., 2002; Peterson et al., 2015), can lead to biases in the estimation of the elastic thickness of the lithosphere.

4 Discussion

4.1 Effect of impurities and liquid water on flexure

Impurities or the presence of liquid water reservoirs can strengthen or weaken the ice shell, which will affect its flexural response to various geological loads. To estimate these effects, we provide two simulated YSEs for Europa. In one case, we added a 500 m layer with a 50% higher yield strength in the middle of the ice shell. The presence of such layer could, for example, result from the existence of a low-porosity or pure-ice region, as the reduction of porosity or impurity content both increase the strength of ice (e.g., John et al., 2018). In a second case, we added a 500 m water ice layer at the ice shell midpoint. For both cases, we compare synthetic flexural profiles to a case without impurities or liquid water (Figure 4).

Figure 4: Effect of ice shell weakening or strengthening on flexural profiles. Yield strength envelopes for the reference, weak, and strong cases (a) and the associated flexural profiles (b).

For the case with impurities strengthening the ice shell, we observe minor effects on the flexural profile. This is due to the additional yield strength only moderately affecting the total bending moment that can be supported by the ice shell (Figure 4b). When including a liquid water pocket with near-zero strength, however, the system becomes mechanically decoupled. Such decoupling drastically reduces the effective elastic thickness and resistance of the plate to yielding (see also Burov, 2010). In that scenario, the flexural profile displays narrower wavelength of flexure with a closer-in flexural bulge. We note that these example models do not require global water ice or impurity layers. Local yield strength anomalies with a wavelength greater than the deformation (i.e. several hundreds of km) would lead to the same flexural profiles. These results imply that if the elastic properties of the ice shell are known, water ice reservoirs and other factors that decouple the flexural response of the ice shell, can be mapped by investigating ice shell flexure. On Earth, elevation change have been typically used to map subglacial lakes in Antartica (Siegert, 2000; Wilson et al., 2025).

4.2 Revisiting flexural profiles at Ariel and Ganymede

Given the potential issues related to purely elastic flexure models, we here revisit two studies that investigated flexure and heat flux at Ariel (Peterson et al., 2015) and Ganymede (Nimmo et al., 2002), using our elastic-plastic model. In Peterson et al. (2015), analyses of steep rift flanks observed in a stereo-photoclinometry digital elevation model constrain the elastic thickness of the lithosphere to be 3.8 to 4.4 km. These values were then used together with the Deborah number formalism to estimate the heat flux, as further described below. Considering a Deborah number of 0.01, three ice rheologies, two end-member strain rate and grain size best-fit heat fluxes at Ariel were found to have been between 28 and 92 mW m\(^{-2}\)  (see also Beddingfield et al., 2022). Using a similar framework, Nimmo et al. (2002) have constrained elastic thicknesses of about 1 km by investigating rift flank uplifts at Ganymede (see also Nimmo & Pappalardo, 2004), and estimated the best-fit heat flux to be \(<\)245 mW m\(^{-2}\).

For both studies, we have digitized the average topography and associated standard deviations profiles of rift flanks considered. Rather than inferring the elastic thickness of the lithosphere from the observed flexural profiles, our approach allows one to directly constrain the yield strength envelope of the ice shell and thus the heat flux for a given rheology. For each case, we compare results using an elastic and elastic-plastic model. Relevant rheological and interior parameters are summarized in Table S1.

4.2.1 The case of Ariel

In Peterson et al. (2015), Ariel’s heat flux is found to have been between 28 and 92 mW m\(^{-2}\) at the time of deformation. For simplicity and comparison purposes, we here focus on the maximum allowed heat flux using our elastic-plastic flexure model. The maximum allowed heat flux is estimated considering the strongest possible ice rheology, which was obtained assuming grain-boundary sliding, a grain-size of 10 cm, and a strain rate of 10\(^{-15}\) s\(^{-1}\) according to the parameter range in Peterson et al. (2015). For consistency, we also assume a thermal conductivity of 5 W m\(^{-1}\) K\(^{-1}\). Similar to the previous section, the brittle part of the ice shell is determined based on parameters found in Beeman et al. (1988) (see also Schulson & Fortt, 2012).

Using these strong rheology parameters, we construct heat flux dependent YSE and estimate elastic-plastic flexure, which is then compared to the broken elastic plate model and observations made in Peterson et al. (2015). For each flexure model, we keep the central deflection \(w_0\) and the reference elevation (the part unaffected by flexure) as free parameters in order to best-fit the observations. We note that additional geometrical parameters could be considered in the fit to observations, such as with the inclusion of a regional slope, and these would affect the inferred heat flux of elastic thickness. However, such considerations are difficult to assess given the limited number of available observations and are beyond the scope of this work. The misfit is determined using a rms considering two cases. In the first case, we only allow models that fit observations within their allowed uncertainty range, which consists of the standard deviation to the mean profile estimated in Peterson et al. (2015) together with an added 150 m uncertainty associated with the vertical resolution of the digital terrain model. In the second case, we retain models with a rms misfit that is less or equal than a total root mean square error encompassing the average standard deviation and elevation uncertainty (both assumed independent).

Misfit plots of heat flux versus \(w_0\) are then constructed for the elastic-plastic model and for the elastic model where we first fit the elastic thickness and then estimate the heat flux using a Deborah number of 0.01. We find that in order to match the observed flexural profiles within uncertainties and with an elastic-plastic model, the heat flux should be less than 32 mW m\(^{-2}\) and with a well-defined best-fit of 25 mW m\(^{-2}\) (Figure 5a,b). For the second approach, the rms misfit is thresholded to 450 m providing a heat flux of at most 45 mW m\(^{-2}\). Weaker rheologies would naturally require lower heat fluxes to match a similar flexural profile. Higher heat fluxes (e.g. 50 mW m\(^{-2}\)) lead to a smaller flexural wavelength that is inconsistent with observations.

When considering an elastic plastic flexure model, where the heat flux is estimated using an assumed Deborah number of 0.01, the allowed heat flux is found to be substantially higher, reaching best-fit values of 92 mW m\(^{-2}\) and allowed values of up to 160 mW m\(^{-2}\) (Figure 5c,d). Thus, when considering elastic-plastic flexure, the required heat flux are substantially lower than found in Peterson et al. (2015) for a similar rheology.

As discussed above, this inconsistency may stem from the assumption of plate elasticity, as the estimated central deflection and plate curvatures are large (\(\sim \)10\(^{-5}\) m\(^{-1}\), Figure S3). Another possibility is related to the Deborah number approach used to estimate the temperature at the base of the elastic lithosphere and the associated heat flux (see Nimmo et al., 2002). The assumed Deborah number of 0.01 represents the materials’ behavior at the crossover depth between the elastic lithosphere and viscous asthenosphere. However, simply translating elastic thickness to elastic depth is an idealized view that does not consider the lithosphere as an elasto-visco-plastic medium. Elasticity within a lithosphere typically occurs within an “elastic-core” located near the midpoint of the plate (e.g., Burov, 2010; McNutt, 1984). Under low curvatures, the entire plate behaves nearly elastically (mechanical thickness equals lithosphere thickness) and the thickness of the lithosphere can be used to estimate the base of the elastic lithosphere. In high curvature regimes, the elastic core is confined to the middle part of the plate, the top of the elastic core does not reach the surface and thus the base of the elastic lithosphere is systematically deeper than its thickness. As a result, the depth of the base of the elastic lithosphere is underestimated using the Deborah number approach and the heat flux is overestimated. This behavior is conceptually observed in the misfit plots, where the Deborah number predicted heat flux increases away from the elastic-plastic fit as the central flexure \(w_0\) (and hence plate curvature) increases (Figures 5 and S3). For small \(w_0\) (\(<\) 3 km) and heat flux (\(<\) 30 mW m\(^{-2}\)), we find curvatures of less than 10\(^{-6}\) m\(^{-1}\) (Figure S3) and both approaches are found to provide comparable fits to observations.

Importantly, we also note that the estimated heat flux is largely affected by unknown parameters controlling the ice shell rheology, such as grain size or strain rate (see also Nimmo et al., 2002; Peterson et al., 2015). These parameters are uncertain by orders of magnitude and can lead to variations in heat flux by a factor of 3 or more (see also Nimmo et al., 2002; Peterson et al., 2015). A full parameter search is beyond the scope of this study, which focuses on comparing elastic-plastic and elastic flexure models.

Figure 5: Rift flank and heat flux at Ariel. Rms misfit contour of heat flux and central deflection \(w_0\) when using our elastic-plastic model to fit the observed flanks at Ariel (a, b). Rms misfit contour using the elastic model and Deborah number approach to estimate the heat flux (c, d). Best-fit flexure profile (black) and rift flank observations with associated uncertainties (blue) digitized from Peterson et al. (2015) (e, f). The black line in a, b, c, d delimits models that fit observations within their 1-sigma bounds. Best-fit parameters are provided in the legend and the red cross indicates the best-fit.

4.2.2 The case of Ganymede

Following a similar approach, we next investigate Ganymede and specifically focus on the rift flanks studied in Nimmo et al. (2002). Before we proceed, we note that the aforementioned study considered two additional geometric parameters to Equation (2), in order to find the best-fit elastic thickness (see also Barnett et al., 2002). However, when only a few observation points are fit and when the observed flexural profile does not extend far enough from the central deflection axis, this updated elastic model formulation leads to the existence of important tradeoffs between the elastic thickness and the geometric parameters. As shown in Figure S2, elastic thicknesses of 0.5 to 1.5 km provide a nearly identical fit to rift flanks at Harpagia Sulcus when using the geometric parameters. Considering a Deborah number of 0.01, grain-boundary sliding, a thermal conductivity of T / 567 and a strain rate of 10\(^{-14}\) s\(^{-1}\), these would imply nearly identical ice shell responses for heat fluxes of 58 to 174 mW m\(^{-2}\). When using Equation (2) and considering the standard deviation to the mean profile estimated in Nimmo et al. (2002), we obtain a best-fit elastic thickness of 0.40 km, with an associated heat flux of 267 mW m\(^{-2}\) (Figure 6a). For comparison purposes and consistency, we consider this latter heat flux to be representative of the deformation at Harpagia Sulcus. Using a similar approach, the heat flux at Arbela Sulcus is found to be 178 mW m\(^{-2}\) (Figures S2 and 6b), which we also retain for comparison to our elastic-plastic model.

Using the same rheological parameters together with our elastic-plastic model, we find that in order to best match the observed rift flank profiles, the best-fit heat flux should be about 163 mW m\(^{-2}\) at Harpagia Sulcus and 117–121 mW m\(^{-2}\) at Arbela Sulcus (Figure 6c,d). Using the first approach to estimate allowable parameters, the heat flux is constrained to be 140 to 205 mW m\(^{-2}\) at Harpagia Sulcus and \(<\)123 mW m\(^{-2}\) at Arbela Sulcus. Using the second and less restrictive approach, models with heat flux up to 300 mW m\(^{-2}\) would be allowed, which results from the large uncertainties in observed topography. As in Nimmo et al. (2002), increasing the strain rate or assuming a weaker rheology would lead to lower required heat flux. As for the case of Ariel, heat fluxes derived using elastic-plastic models are systematically lower than estimated using the elastic model and Deborah number approach, though we note that solutions overlap substantially given the large uncertainties in the observed profiles.

Figure 6: Rift flank and heat flux at Ganymede. Rms misfit contour of heat flux and central deflection \(w_0\) when using our elastic-plastic model to fit the observed flanks at Ganymede (a, b). Rms misfit contour using the elastic model and Deborah number approach to estimate the heat flux (c,d). Best-fit flexural profile (red), best-fit profiles that fit observations within their allowed uncertainties (black), and rift flank observations with associated uncertainties (blue) digitized from Nimmo et al. (2002) for Harpagia Sulcus (e) and Arbela Sulcus (f). The black line in a, b, c, d delimit models that fit observations within their 1-sigma bounds with the red cross showing the associated best-fit. The black dot shows the best-fit for any model. Best-fit parameters are provided in the legend.

4.3 On the static assumption

Our ice-shell flexure models make the important assumption that the observed deformation is static and/or at steady state, which is only valid if the time elapsed since the deformation started is greater than the time required for viscous adjustments (e.g., Peltier, 1974). Icy moons are dynamic worlds (e.g., Dombard & McKinnon, 2006; McKinnon, 1998; Van Hoolst et al., 2024) and it is likely that some of the observed surface deformations are still ongoing. However, considering the ice-shell transient viscoelastic response would require knowledge on the age of the observed deformation features, which is currently poorly constrained (e.g., Čadek et al., 2017). Nevertheless, it is important to note that modeling the ice-shell transient response with a static model in which deformations are instantaneous would generally lead to an overestimation of the ice-shell yield strength (and underestimation of the heat flux). Indeed, if the deformation is still ongoing, only a fraction of the final deformation would be observed and this would make a static model misinterpret the shell as being stronger and colder than it is. Any subsequent viscoelastic relaxation processes unrelated to the loading event would also tend to alter and reduce the observed deformation wavelengths and lead to an overestimation of the ice-shell yield strength. We note that such problems related to estimates of paleo elastic flexure and heat flux are common in geophysical modeling of the terrestrial planets (e.g., Broquet & Wieczorek, 2019; Smrekar et al., 2022).

Importantly, as the predicted heat fluxes for Ganymede are far above the estimated radiogenic heat fluxes of 40 mW m\(^{-2}\) (e.g., Bland et al., 2009), it is likely that the observed deformations have been recorded during a time discrete major tidal heating event. Such deformations are thus likely to have been frozen in the ice shell and be preserved to present-day.

5 Conclusion

Investigations of surface flexure provide critical insights to the interior structure and thermal state of planetary bodies. In most previous works looking into icy worlds, observed flexural signatures at rift zones or furrows were fitted using an analytical broken elastic plate model in order to estimate the elastic thickness of the lithosphere and associated heat flux (e.g., Nimmo et al., 2002; Peterson et al., 2015). In this work, we discuss a more adequate estimation of flexure that considers ice shell plasticity. Rather than estimating the elastic thickness of the lithosphere and associated heat flux, the presented elastic-plastic model allows to directly consider the effect of ice-rheology and heat flux on the flexural response of the ice shell. Such model is not computationally expensive to run, does not rely on the elastic thickness mathematical construct, and should thus be preferred over the classical broken elastic plate model. Nevertheless, we find that elastic and elastic-plastic solutions can agree well in low plate curvature geologic settings, such that elastic models remain adequate below a specific plate rheology dependent bending-moment saturation threshold. We also note that the model does not account for the possibility that observed deformations are undergoing viscoelastic relaxation, which would affect the inferred heat flux.

Comparing elastic-plastic to broken elastic plate models, we show that the maximum curvature of the synthetic model should be used when relating an elastic thickness estimate to a yield strength envelope through the equating bending moment approach. We further identify that in regions where the plate curvature exceeds a critical threshold, the plate yields making its flexural response strongly depart from analytical broken elastic plate model predictions. This effect is prominently expressed before the first zero-crossing point and in the flexural bulge region, where the elastic model predicts unrealistic surface oscillations. Our analysis implies that elastic models are adequately suited to estimate the elastic thickness of the lithosphere and associated heat flux, as long as the observed deformed surface displays curvatures that are below a saturation point. At such saturation level, any increase in plate curvature should not be associated with an increase in bending moment, a behavior that is not considered by elastic flexure theory. The saturation point of the ice shell is a function of its yield strength envelope and can be determined numerically. However, as the ice shell curvature in elastic models depend on both the elastic thickness and central deflection value, it is not possible to provide an elastic thickness threshold below which purely elastic models become unreliable.

Local strengthening on the ice shell due, for example, to a reduction of the fraction of impurities was found to have little effect on the predicted plate flexure. On the other hand, the presence of a near-zero yield strength zone in the ice shell can decouple its mechanical response to loading, reduce its effective elastic thickness and affect its flexural response. Thus the presence of local weak zones, such as liquid water pockets, can potentially be detected in elevation data assuming the elastic thickness or ice yield strength is adequately known, as has been done on Earth (Siegert, 2000; Wilson et al., 2025). Comparisons of flexural profiles in various locations and where ice rheology should remain constant could help detecting such ice shell layering and could be used in conjunction with ice penetrating radars (e.g., Blankenship et al., 2024; Bruzzone et al., 2013).

We have applied our elastic-plastic flexure model to rift flanks at Ganymede and Ariel, focusing specifically on regions investigated in the works of Nimmo et al. (2002) and Peterson et al. (2015). We find the heat flux, in these high curvature geologic settings \(>\)10\(^{-6}\) m\(^{-1}\), to have been overestimated. The overestimation results from both the assumption of pure elasticity and the inadequate estimation of the heat flux using a Deborah number approach in high plate curvature settings. Our results are generally valid for any analysis of flexural profiles, and have implications for our understanding of the geologic history of icy moons.

Open science statements

Author contributions A.B. is responsible for the conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing (original draft), and writing (review & editing).

Code availability Python-based codes to compute elastic-plastic flexure, yield strength envelope and heat fluxes can be found at Broquet (2022) (https://doi.org/10.5281/zenodo.4973893) and Broquet (2026) (https://doi.org/10.5281/zenodo.17047057).

Data availability Europa’s example yield strength envelope together with digitized elevation data for Ganymede and Ariel can be found at Broquet (2026) (https://doi.org/10.5281/zenodo.17047057).

Funding A.B. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for his Emmy Noether junior research group – Project-ID 559317179.

Competing interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

Barnett, D. N., Nimmo, F., & McKenzie, D. (2002). Flexure of venusian lithosphere measured from residual topography and gravity. Journal of Geophysical Research: Planets, 107(E2), 2-1-2–21. https://doi.org/https://doi.org/10.1029/2000JE001398

Beddingfield, C. B., Cartwright, R. J., Leonard, E., Nordheim, T., & Scipioni, F. (2022). Ariel’s elastic thicknesses and heat fluxes. The Planetary Science Journal, 3(5), 106. https://doi.org/10.3847/psj/ac63d1

Beeman, M., Durham, W. B., & Kirby, S. H. (1988). Friction of ice. Journal of Geophysical Research: Solid Earth, 93(B7), 7625–7633. https://doi.org/10.1029/jb093ib07p07625

Berne, A., Simons, M., Keane, J. T., Leonard, E. J., & Park, R. S. (2024). Jet activity on enceladus linked to tidally driven strike-slip motion along tiger stripes. Nature Geoscience, 17(5), 385–391. https://doi.org/10.1038/s41561-024-01418-0

Billings, S. E., & Kattenhorn, S. A. (2005). The great thickness debate: Ice shell thickness models for europa and comparisons with estimates based on flexure at ridges. Icarus, 177(2), 397–412. https://doi.org/10.1016/j.icarus.2005.03.013

Bland, M. T., Showman, A. P., & Tobie, G. (2009). The orbital–thermal evolution and global expansion of ganymede. Icarus, 200(1), 207–221. https://doi.org/10.1016/j.icarus.2008.11.016

Blankenship, D. D., Moussessian, A., Chapin, E., Young, D. A., Wesley Patterson, G., Plaut, J. J., Freedman, A. P., Schroeder, D. M., Grima, C., Steinbrugge, G., Soderlund, K. M., Ray, T., Richter, T. G., Jones-Wilson, L., Wolfenbarger, N. S., Scanlan, K. M., Gerekos, C., Chan, K., Seker, I., … Urban, T. J. (2024). Radar for europa assessment and sounding: Ocean to near-surface (reason). Space Science Reviews, 220(5). https://doi.org/10.1007/s11214-024-01072-3

Broquet, A. (2022). AB-Ares/Te_HF_Conversion: 0.2.3. https://doi.org/10.5281/zenodo.4973893

Broquet, A. (2026). Data and code for elastic-plastic flexure on icy moons. https://doi.org/10.5281/zenodo.17047057

Broquet, A., & Wieczorek, M. A. (2019). The gravitational signature of martian volcanoes. Journal of Geophysical Research: Planets, 124(8), 2054–2086. https://doi.org/https://doi.org/10.1029/2019JE005959

Bruzzone, L., Plaut, J. J., Alberti, G., Blankenship, D. D., Bovolo, F., Campbell, B. A., Ferro, A., Gim, Y., Kofman, W., Komatsu, G., McKinnon, W. B., Mitri, G., Orosei, R., Patterson, G. W., Plettemeier, D., & Seu, R. (2013). Rime: Radar for icy moon exploration. 2013 IEEE International Geoscience and Remote Sensing Symposium - IGARSS, 3907–3910. https://doi.org/10.1109/igarss.2013.6723686

Buffo, J. J., Schmidt, B. E., Huber, C., & Walker, C. C. (2020). Entrainment and dynamics of ocean‐derived impurities within europa’s ice shell. Journal of Geophysical Research: Planets, 125(10). https://doi.org/10.1029/2020je006394

Burov, E. B. (2010). The equivalent elastic thickness (te), seismicity and the long-term rheology of continental lithosphere: Time to burn-out “crème brûlée”? Tectonophysics, 484(1–4), 4–26. https://doi.org/10.1016/j.tecto.2009.06.013

Burov, E. B., & Diament, M. (1995). The effective elastic thickness (te) of continental lithosphere: What does it really mean? Journal of Geophysical Research: Solid Earth, 100(B3), 3905–3927. https://doi.org/10.1029/94jb02770

Čadek, O., Běhounková, M., Tobie, G., & Choblet, G. (2017). Viscoelastic relaxation of enceladus’s ice shell. Icarus, 291, 31–35. https://doi.org/10.1016/j.icarus.2017.03.011

Damptz, A. L., & Dombard, A. J. (2011). Time-dependent flexure of the lithospheres on the icy satellites of jupiter and saturn. Icarus, 216(1), 86–88. https://doi.org/10.1016/j.icarus.2011.07.011

Dombard, A. J., & McKinnon, W. B. (2006). Folding of europa’s icy lithosphere: An analysis of viscous-plastic buckling and subsequent topographic relaxation. Journal of Structural Geology, 28(12), 2259–2269. https://doi.org/10.1016/j.jsg.2005.12.003

Enya, K., Kobayashi, M., Kimura, J., Araki, H., Namiki, N., Noda, H., Kashima, S., Oshigami, S., Ishibashi, K., Yamawaki, T., Tohara, K., Saito, Y., Ozaki, M., Mizuno, T., Kamata, S., Matsumoto, K., Sasaki, S., Kuramoto, K., Sato, Y., … Hussmann, H. (2022). The ganymede laser altimeter (gala) for the jupiter icy moons explorer (juice): Mission, science, and instrumentation of its receiver modules. Advances in Space Research, 69(5), 2283–2304. https://doi.org/10.1016/j.asr.2021.11.036

Fung, Y. C. (1965). Foundations of solid mechanics. Prentice-Hall, Englewood Cliffs, NJ.

Giese, B., Wagner, R., Hussmann, H., Neukum, G., Perry, J., Helfenstein, P., & Thomas, P. C. (2008). Enceladus: An estimate of heat flux and lithospheric thickness from flexurally supported topography. Geophysical Research Letters, 35(24). https://doi.org/10.1029/2008gl036149

Goetze, C., & Evans, B. (1979). Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics. Geophysical Journal International, 59(3), 463–478. https://doi.org/10.1111/j.1365-246x.1979.tb02567.x

Grasset, O., Dougherty, M., Coustenis, A., Bunce, E., Erd, C., Titov, D., Blanc, M., Coates, A., Drossart, P., Fletcher, L., Hussmann, H., Jaumann, R., Krupp, N., Lebreton, J.-P., Prieto-Ballesteros, O., Tortora, P., Tosi, F., & Van Hoolst, T. (2013). Jupiter icy moons explorer (juice): An esa mission to orbit ganymede and to characterise the jupiter system. Planetary and Space Science, 78, 1–21. https://doi.org/10.1016/j.pss.2012.12.002

Grott, M., Hauber, E., Werner, S. C., Kronberg, P., & Neukum, G. (2005). High heat flux on ancient mars: Evidence from rift flank uplift at coracis fossae. Geophysical Research Letters, 32(21). https://doi.org/10.1029/2005gl023894

Hammond, N. P., Phillips, C. B., Nimmo, F., & Kattenhorn, S. A. (2013). Flexure on dione: Investigating subsurface structure and thermal history. Icarus, 223(1), 418–422. https://doi.org/10.1016/j.icarus.2012.12.021

Hurford, T. A., Beyer, R. A., Schmidt, B., Preblich, B., Sarid, A. R., & Greenberg, R. (2005). Flexure of europa’s lithosphere due to ridge-loading. Icarus, 177(2), 380–396. https://doi.org/10.1016/j.icarus.2005.06.019

John, M., Suominen, M., Sormunen, O.-V., Hasan, M., Kurvinen, E., Kujala, P., Mikkola, A., & Louhi-Kultanen, M. (2018). Purity and mechanical strength of naturally frozen ice in wastewater basins. Water Research, 145, 418–428. https://doi.org/10.1016/j.watres.2018.08.063

Mancktelow, N. S. (1999). Finite-element modelling of single-layer folding in elasto-viscous materials: The effect of initial perturbation geometry. Journal of Structural Geology, 21(2), 161–177. https://doi.org/10.1016/s0191-8141(98)00102-3

McKinnon, W. B. (1998). Geodynamics of icy satellites. In Solar system ices (pp. 525–550). Springer Netherlands. https://doi.org/10.1007/978-94-011-5252-5˙22

McNutt, M. K. (1984). Lithospheric flexure and thermal anomalies. Journal of Geophysical Research: Solid Earth, 89(B13), 11180–11194. https://doi.org/10.1029/jb089ib13p11180

McNutt, M. K., & Menard, H. W. (1982). Constraints on yield strength in the oceanic lithosphere derived from observations of flexure. Geophysical Journal International, 71(2), 363–394. https://doi.org/10.1111/j.1365-246x.1982.tb05994.x

Mueller, S., & Phillips, R. J. (1995). On the reliability of lithospheric constraints derived from models of outer-rise flexure. Geophysical Journal International, 123(3), 887–902. https://doi.org/10.1111/j.1365-246x.1995.tb06896.x

Nimmo, F. (2018). Icy satellites: Interior structure, dynamics, and evolution. https://doi.org/10.1093/acrefore/9780190647926.013.29

Nimmo, F., & Pappalardo, R. T. (2004). Furrow flexure and ancient heat flux on ganymede. Geophysical Research Letters, 31(19). https://doi.org/10.1029/2004gl020763

Nimmo, F., Pappalardo, R. T., & Giese, B. (2002). Effective elastic thickness and heat flux estimates on ganymede. Geophysical Research Letters, 29(7). https://doi.org/10.1029/2001gl013976

Park, R. S., Mastrodemos, N., Jacobson, R. A., Berne, A., Vaughan, A. T., Hemingway, D. J., Leonard, E. J., Castillo‐Rogez, J. C., Cockell, C. S., Keane, J. T., Konopliv, A. S., Nimmo, F., Riedel, J. E., Simons, M., & Vance, S. (2024). The global shape, gravity field, and libration of enceladus. Journal of Geophysical Research: Planets, 129(1). https://doi.org/10.1029/2023je008054

Peltier, W. R. (1974). The impulse response of a maxwell earth. Reviews of Geophysics, 12(4), 649–669. https://doi.org/10.1029/rg012i004p00649

Peterson, G., Nimmo, F., & Schenk, P. M. (2015). Elastic thickness and heat flux estimates for the uranian satellite ariel. Icarus, 250, 116–122. https://doi.org/10.1016/j.icarus.2014.11.007

Phillips, R. J. (1990). Convection‐driven tectonics on venus. Journal of Geophysical Research: Solid Earth, 95(B2), 1301–1316. https://doi.org/10.1029/jb095ib02p01301

Schenk, P. M., & McKinnon, W. B. (2024). New global topography of enceladus: Hypsometry, basins, spherical harmonics, shell thickness, and true polar wander revisited. Icarus, 408, 115827. https://doi.org/10.1016/j.icarus.2023.115827

Schulson, E. M., & Fortt, A. L. (2012). Friction of ice on ice. Journal of Geophysical Research: Solid Earth, 117(B12). https://doi.org/10.1029/2012jb009219

Siegert, M. J. (2000). Antarctic subglacial lakes. Earth-Science Reviews, 50(1–2), 29–50. https://doi.org/10.1016/s0012-8252(99)00068-9

Smrekar, S. E., Ostberg, C., & O’Rourke, J. G. (2022). Earth-like lithospheric thickness and heat flow on venus consistent with active rifting. Nature Geoscience, 16(1), 13–18. https://doi.org/10.1038/s41561-022-01068-0

Souček, O., Běhounková, M., Schroeder, D. M., Wolfenbarger, N. S., Kalousová, K., Steinbrugge, G., & Soderlund, K. M. (2023). Radar attenuation in enceladus’ ice shell: Obstacles and opportunities for constraining shell thickness, chemistry, and thermal structure. Journal of Geophysical Research: Planets, 128(2). https://doi.org/10.1029/2022je007626

Spitale, J. N., Tigges, M. D., Berne, A., Rhoden, A., Hurford, T. A., & Webster, K. D. (2025). Curtain-based maps of eruptive activity in enceladus’s south-polar terrain at 15 cassini epochs. The Planetary Science Journal, 6(3), 67. https://doi.org/10.3847/psj/adb7d7

Turcotte, D. L., & Schubert, G. (2002). Geodynamics. Cambridge University Press. https://doi.org/10.1017/cbo9780511807442

Van Hoolst, T., Tobie, G., Vallat, C., Altobelli, N., Bruzzone, L., Cao, H., Dirkx, D., Genova, A., Hussmann, H., Iess, L., Kimura, J., Khurana, K., Lucchetti, A., Mitri, G., Moore, W., Saur, J., Stark, A., Vorburger, A., Wieczorek, M., … Witasse, O. (2024). Geophysical characterization of the interiors of ganymede, callisto and europa by esa’s jupiter icy moons explorer. Space Science Reviews, 220(5). https://doi.org/10.1007/s11214-024-01085-y

Vance, S. D., Craft, K. L., Shock, E., Schmidt, B. E., Lunine, J., Hand, K. P., McKinnon, W. B., Spiers, E. M., Chivers, C., Lawrence, J. D., Wolfenbarger, N., Leonard, E. J., Robinson, K. J., Styczinski, M. J., Persaud, D. M., Steinbrugge, G., Zolotov, M. Y., Quick, L. C., Scully, J. E. C., … Elder, C. M. (2023). Investigating europa’s habitability with the europa clipper. Space Science Reviews, 219(8). https://doi.org/10.1007/s11214-023-01025-2

Wilson, S. F., Hogg, A. E., Rigby, R., Gourmelen, N., Nias, I., & Slater, T. (2025). Detection of 85 new active subglacial lakes in antarctica from a decade of cryosat-2 data. Nature Communications, 16(1). https://doi.org/10.1038/s41467-025-63773-9

Woltring, H. J. (1986). A fortran package for generalized, cross-validatory spline smoothing and differentiation. Advances in Engineering Software (1978), 8(2), 104–113. https://doi.org/10.1016/0141-1195(86)90098-7