Introduction
Grounding zones are transition regions where grounded ice sheets become floating ice shelves and ice, ocean and the solid Earth all meet. The dynamics of the grounding zone are among the most sensitive indicators of change in the Antarctic (Pattyn, Reference Pattyn2017; Scambos and others, Reference Scambos2017; Gudmundsson and others, Reference Gudmundsson, Paolo, Adusumilli and Fricker2019; Seroussi and others, Reference Seroussi2020). The satellite era has revealed widespread, rapid grounding line retreat (Joughin and others, Reference Joughin, Smith and Medley2014; Rignot and others, Reference Rignot, Mouginot, Morlighem, Seroussi and Scheuchl2014), especially in West Antarctica, where ice shelves are increasingly melted from below by warm, salty Circumpolar Deep Water (Jourdain and others, Reference Jourdain2017; Nakayama and others, Reference Nakayama, Menemenlis, Zhang, Schodlok and Rignot2018). Basal melt rate in the grounding zone is one of the largest sources of uncertainty in future projections of the rate and amount of future sea-level rise (Pritchard and others, Reference Pritchard, Ligtenberg, Fricker, Vaughan, van den Broeke and Padman2012; Adusumilli and others, Reference Adusumilli, Fricker, Medley, Padman and Siegfried2020; Seroussi and others, Reference Seroussi2020; Hill and others, Reference Hill, Rosier, Gudmundsson and Collins2021).
Basal melt rate 
$m$ can be calculated by conservation of mass, such that:
\begin{equation}
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\frac{\partial h}{\partial t} = (a - m) - \nabla \cdot (h\underline{u})
\end{equation}where 
$a$ is the surface mass balance and 
$\underline{u}$ is the ice velocity, provided that the local ice thickness 
$h$ is well-known. Ice shelf thickness is typically calculated by a freeboard approach, where the ice height above the ocean surface is measured and the total ice thickness is calculated as the freeboard height times 
$\rho_{sea}/(\rho_{sea} - \rho_{ice})$, with a correction for firn densification, for example, Smith and others (Reference Smith2020); Chartrand and Howat (Reference Chartrand and Howat2023). This works well away from the grounding line, where ice is in hydrostatic equilibrium and floats up and down on ocean tides. However, near the grounding line, the ice shelf is mechanically coupled to the upstream ice sheet and flexes, rather than floats (Fig. 1).
Figure 1. We model the flexure zone as an elastic beam under tidal forcing with effective Young’s modulus E*. The far-field sea level is 
$A_0(t)$. At 
$A_0$ = 0, 
$w(x)$ = 0, and the neutral surface of our model beam rests on the 
$x$-axis. As the tide changes, the upward force on the beam is the hydrostatic pressure proportional to the difference between 
$A_0$ and 
$w(x)$. We can observe this flexure by differencing repeat track ice elevation measurements from all available ICESat-2 cycles. We allow the ice thickness in the grounding zone to vary as 
$h(x)$ and model the resultant combined flexure of such a beam with an effective Young’s modulus E*. Surface and basal crevasses may be present. Tidal mixing takes place at the ice–ocean interface, and the grounding line may move back and forth. ICESat-2 image adapted from National Aeronautics and Space Administration (2025). Schematic illustration of a survey aircraft (aircraft silhouette, CC0 from Wikimedia Commons).
The flexure of an elastic beam under a load depends on its flexural rigidity, proportional to the product of its thickness cubed and Young’s modulus 
$E$ (Turcotte and Schubert, Reference Turcotte and Schubert2014). Vaughan (Reference Vaughan1995) modeled the tidal flexure of ice as an elastic beam of constant thickness and measured vertical tidal deflection on the Rutford Ice Stream with repeated surveys of the same transect at different points in the tidal cycle with a Global Positioning System (GPS)-equipped snowmobile. By inverting the observable flexure for the properties of the beam, Vaughan found that a value for 
$E$ of 0.88 
$\pm$ 0.35 GPa adequately minimized the misfit between modeled and observed flexure using a constant ice thickness 
$h$ from contemporaneous radar sounding measurements, as only one of 
$h$ or 
$E$ can be inferred given the other.
This result is somewhat surprising, as the value of the Young’s modulus of ice is commonly taken to be 9 GPa, from laboratory measurements: an order of magnitude greater (Hooke, Reference Hooke2005). The Young’s modulus of a material is an elastic parameter that reflects the ratio of an applied uniaxial normal stress to instantaneous normal strain and has its physical basis in the electromagnetic forces between molecules in its crystal lattice. It is measured in the laboratory by mechanical stress–strain curve testing, or by wave speed measurement through the lattice, or by interferometric phonon scattering measurements (Gammon and others, Reference Gammon, Kiefte, Clouter and Denner1983; Schulson and Duval, Reference Schulson and Duval2011; Rathmann and others, Reference Rathmann2022). As ice is a viscoelastic material, there is a dependence of E on the rate of the applied stress due to creep and anelastic effects happening over time, as well as on temperature and presence of water or inclusions (Sinha, Reference Sinha1978; Reference Sinha1989). Sinha (Reference Sinha1978) uses the term effective Young’s modulus to distinguish between the ‘true’ instantaneous elastic modulus at high frequencies and the relationship between stress and strain of in situ ice. Gold (Reference Gold1977) compiles a number of laboratory experiments and finds an approximately linear relationship between Young’s modulus and the frequency of applied stress, leading to E 
$\approx$ 3 GPa at 
$10^{-3}$ Hz. The tidal frequency is about 
$10^{-5}$ Hz, which would give E 
$\approx$ 1 GPa if the linear relationship continues, as noted by Vaughan (Reference Vaughan1995) and a factor of ten from the textbook value of 9 GPa.
Numerous approaches have been employed to resolve the apparent discrepancy between the laboratory values of E = 9 GPa and the much smaller results from tidal flexure models. Reeh and others (Reference Reeh, Christensen, Mayer and Olesen2003) develop a linear viscoelastic model similar to a Burgers model with four elastic parameters and use flexure data on Nioghalvfjerdsfjorden glacier from tiltmeters and GPS, with airborne ice-penetrating radar thickness, and temperature measurements from a borehole to determine the local ice viscosity. Detailed information on the tidal phase is also needed, which is possible to acquire with tiltmeters (e.g., Smith, Reference Smith1991) but not with satellite laser altimetry. The launch of ICESat (Ice, Cloud and land Elevation Satellite) in 2003, followed by ICESat-2 in 2018, made it possible to difference repeat track measurements of ice surface height at different points in the tidal cycle to observe tidal flexure near the grounding line on many ice shelves around Antarctica (Fricker and Padman, Reference Fricker and Padman2006). Bindschadler and others (Reference Bindschadler2011) mapped grounding zones around Antarctica with ICESat for the International Polar Year 2007–09 and applied the model of Vaughan (Reference Vaughan1995) to compare the observed and expected relationship between flexure zone length and ice thickness, but found considerable scatter. Sayag and Worster (Reference Sayag and Worster2013) modeled the flexure as an elastic beam of constant thickness with varied bedrock and grounding line conditions combined with ICESat altimetry in two places on the Ronne-Filchner Ice Shelf. They find a best fitting value for E of 1.8 GPa in a ‘stiff-fixed’ model, where the grounding line is fixed and ice rests on a stiff frozen bed, but a best fitting E of 9.33 GPa with a ‘soft-free’ model, where the grounding line is allowed to move and ice rests on an elastic bed having an estimated spring constant of 1 MPa/m. Marsh and others (Reference Marsh, Rack, Golledge, Lawson and Floricioiu2014) used differential interferometric synthetic aperture radar (InSAR) to observe tidal flexure on Beardmore Glacier. They allow the ice thickness to vary and, using a value for E of 1.4 GPa, are able to invert for local ice thickness using COMSOL, a commercial multi-physics finite-element software package.
A number of studies have incorporated viscoelastic effects and compared these to elastic models. Wild and others (Reference Wild, Marsh and Rack2018) compare a viscoelastic and an elastic model, with an elastic bed of spring constant 5 MPa/m, and find that the viscoelastic model is very sensitive to the tidal model used to determine the phase of the tide. Rosier and others (Reference Rosier, Marsh, Rack, Gudmundsson, Wild and Ryan2017) compare a simple elastic beam model to full-Stokes viscoelastic models with inclusion of basal crevasses and density-dependent ice stiffness and find that the elastic solution model can produce similarly excellent fits to data. They also find that a factor of two of variation, ranging from about 2 to 4 GPa, in estimated effective Young’s modulus exists even when local ice thickness is well constrained.
The in situ environment of ice near the grounding line is complex and influenced by physical forces on intersecting length and timescales (Fig. 1). Ocean tides cycle around Antarctica every 12–24 hours (Padman and others, Reference Padman, Siegfried and Fricker2018). Tides pump seawater in and out of the flexure zone, acting as a forcing to the hydrodynamics of a poroelastic bed with sediment on a spectrum from fully frozen to fully deformable (Warburton and others, Reference Warburton, Hewitt and Neufeld2020; Kazmierczak and others, Reference Kazmierczak, Gregov, Coulon and Pattyn2024). There is tidal mixing of incoming and outgoing water masses: warm, salty Circumpolar Deep Water increasingly encroaches on the Antarctic continental shelf and melts ice from below, releasing cold glacial meltwater, which can refreeze onto the ice shelf (Nakayama and others, Reference Nakayama, Menemenlis, Zhang, Schodlok and Rignot2018). There are surface and basal crevasses and channels that direct the flow of water (Alley and others, Reference Alley, Scambos and Alley2022). Recent direct underwater observations of the grounding zone have begun to show a complex system where these and other processes interact on different length and timescales to produce angular terraced forms on the undersides of ice shelves (Dutrieux and others, Reference Dutrieux2014), as well as swooping, curved regions where convective and oceanographic forcings act differently on the flexure zone system (Schmidt and others, Reference Schmidt2023).
The observable surface flexure of ice near the grounding line reflects the combined effects of these competing processes, as well as viscous and anelastic processes that affect the stress–strain relationship in real time. For this reason, we will refer only to an effective Young’s modulus E*, to reflect the fact that the relationship between stress and strain that can be inferred from tidal flexure is a bulk parameter not identical to the Young’s modulus of ice as can be measured in the lab, and may indeed vary considerably in space.
Here, we use satellite observations of the surface flexure of ice along with a simple physical model of ice in the flexure zone as an elastic beam of varying thickness, motivated by the need for better estimates of ice thickness and thickness gradient near the grounding line, to constrain the effective Young’s modulus of ice in three places on the Ross Ice Shelf (RIS). We infer a single effective elastic parameter, the effective Young’s modulus E*, to parameterize the observable flexure of the combined ice–bed–ocean system near the grounding line, and argue that we cannot readily distinguish flexure of the ice from flexure of the bed. This allows much of the uncertainty inherent to the problem to be wrapped into one effective parameter whose associated uncertainty can be quantified. We use airborne ice-penetrating radar measurements of ice thickness 
$h(x)$ on the RIS (Das and others, Reference Das2020) to constrain the effective rheology of ice in the flexure zone by minimizing the misfit between modeled flexure and observed ice flexure with laser altimetry data from ICESat-2.
Methods
We employ an elastic beam bending model constrained by repeat track flexure data from ICESat-2 to estimate the flexural rigidity of ice in the flexure zone at three sites on the RIS. These sites are chosen to be far from confining topography, with ICESat-2 ground tracks roughly perpendicular to the grounding line and close to ice thickness measurements, in order to narrowly apply the linear elastic approximation. Using ice-penetrating radar thickness data from the ROSETTA-Ice airborne geophysical survey (Das and others, Reference Das2020), we then calculate the effective Young’s modulus of ice in the flexure zone.
Model
We model ice in the flexure zone as an elastic beam of varying thickness, under small deformations due to tidal forcing, such that the linear elastic approximation holds and the resultant vertical deflection 
$w(x)$ is described by the Euler-Bernoulli beam bending equation, after Holdsworth (Reference Holdsworth1969) and Robin (Reference Robin1958):
\begin{equation}
\frac{d^2}{dx^2}\Big[D(x)\frac{d^2 w}{dx^2}\Big] = \rho_{sea} g \big[A_0 - w(x)\big]
\end{equation}where 
$A_0$ is the far-field sea level, 
$\rho_{sea}$ is the mass density of seawater, 
$g$ is gravitational acceleration and 
$D$ is the spatially variable flexural rigidity of the beam:
\begin{equation}
D(x) = \frac{Eh^3}{12(1-\nu^2)}
\end{equation}where 
$\nu$ is the Poisson’s ratio of ice, here taken to be 
$0.3$ (Schulson, Reference Schulson1999).
To solve, we discretize (2) by central differences after Jacquot and Dewey (Reference Jacquot and Dewey2001) (see Supplementary Material S1). Given full information about 
$D$, that is, 
$h$ and 
$E$, the flexure of a beam under an applied forcing can be readily calculated. This is called the forward problem. The inverse problem seeks to infer properties of the beam given the flexure 
$w$. However, we will only ever be able to invert or optimize for the flexural rigidity 
$D$ of the beam, that is, only one of the Young’s modulus or the thickness, assuming the other is known. This problem is also ill-posed, making it sensitive to noise and resulting in solutions that may oscillate about the correct one (Lucchinetti and Stüssi, Reference Lucchinetti and Stüssi2002).
Solving (2) requires four boundary conditions. The model results are sensitive to these boundary conditions and where they are applied. Here, we select our problem area as the flexure zone: everywhere flexure is observed to occur on tidal timescales. The landward boundary is where vertical tidal deflection vanishes. The seaward boundary is where the vertical tidal deflection reaches the constant value of 
$A_0$, the far-field sea level, which varies at each repeat track measurement. That is, 
$\partial w / \partial x = 0$ and 
$w = 0$ at 
$x = 0$, and 
$\partial w / \partial x = 0 $ and 
$ w = A_0$ at 
$ x = L$, where 
$L$ is the length of the flexure zone. This is the same set of boundary conditions used by Vaughan (Reference Vaughan1995). The length of the flexure zone and 
$A_0$ are inferred from ICESat-2 data, as described in the ICESat-2 grounding zone deflection data section.
Then, using spatially varying ice thickness from deep ice radar (DICE) (see below), we use the MATLAB patternsearch solver (The MathWorks Inc., 2022), a global optimization algorithm, to solve (2) many (e.g., hundreds or thousands) of times while varying E* to find the value of E* that minimizes the square of the misfit between the model and data. In order to manage the sensitivity of the problem to the location of the boundaries, we repeat the optimization with the landward boundary at a range of discrete points up to 1 km upstream and downstream of its initial estimated position. This is akin to a free grounding line position, though we note that the landward extent of tidal flexure may well be upstream of the grounding line.
Our method includes a grid search (see subsection below) that sweeps across the parameter space of E* and the position of the landward boundary, seeking a minimum in the misfit space. For each beam and track at each site, we infer a value for E*. As the problem is ill-posed, each solution to the inverse problem is checked against several criteria before including it as a result. These are discussed in detail in Supplementary Material S1, along with all included and excluded results from the inversion.
Variable ice shelf thickness data
Here we use ice thickness data from the ROSETTA-Ice airborne survey from 2015 to 2017 (Das and others, Reference Das2020), which uses deep ice radar (DICE) with a 188 MHz center frequency, 60 MHz bandwidth and 1.4 m range resolution to produce a dataset of ice thickness on the RIS (Fig. 2). To find spatially varying ice thickness from DICE at points that are coincident with ICESat-2 ground tracks in our selected regions of interest, we use all the DICE datapoints as a scattered interpolant in order to linearly interpolate between them.
Figure 2. The RIS and its location on the Antarctic continent (inset). The three regions studied here are highlighted: Siple Coast (SC), MacAyeal Ice Stream (MIS) and Marie Byrd Land (MBL). Study sites were chosen for their proximity to DICE measurements, distance from confining topography and flexure data consistent with the boundary conditions we apply here. Horizontal and vertical lines show flight paths from the ROSETTA-Ice airborne ice-penetrating radar campaign (Das and others, Reference Das2020). Colors in panel (a) show ice surface velocity (Mouginot and others, Reference Mouginot, Scheuchl and Rignot2017). Panel (b) shows the average inferred effective Young’s modulus along each beam pair ground track at each site. Sections of the ICESat-2 ground tracks used for modeling are shown, dotted. Coordinates in (a) are provided at selected points for spatial reference. A 10 km scale bar is shown in (b).
Interpreting radargrams in the grounding zone is more difficult than it is near the calving front because seawater intrusion upstream of the grounding line and grounding zone hydrodynamics can complicate delineation of the bottom surface of the ice (MacGregor and others, Reference MacGregor, Anandakrishnan, Catania and Winebrenner2011). While crossover analysis in Das and others (Reference Das2020) was robust over the 3 years of data collection and consistent to 2 m, the uncertainty in ice thickness in the grounding zone is likely greater, though these are among the most detailed direct inferences of ice thickness that presently exist. Continent-wide data products of ice shelf thickness typically calculate ice thickness by assuming that ice downstream of a fixed grounding line is free-floating and in hydrostatic equilibrium, then make a correction based on any existing nearby radar ice thickness measurements and interpolate (Morlighem and others, Reference Morlighem2020). This is an excellent assumption, far from the grounding line that breaks down in the grounding zone, where ice is mechanically coupled to the grounded ice sheet. In places, sub-ice-shelf bathymetry measurements are sparse, and uncertainty in ice shelf thickness may exceed 500 m (Morlighem and others, Reference Morlighem2020).
Boundary value grid search
The location of the landward limit of tidal flexure, 
$x_F$, can migrate several kilometers on tidal timescales, as is seen in both observations and models, for example, Robel and others (Reference Robel, Wilson and Seroussi2022); Freer and others (Reference Freer, Marsh, Hogg, Fricker and Padman2023); Stubblefield and others (Reference Stubblefield, Meyer, Siegfried, Sauthoff and Spiegelman2023), and is often identified as Point F in studies of tidal flexure, for example, Fricker and Padman (Reference Fricker and Padman2006); Brunt and others (Reference Brunt, Fricker, Padman, Scambos and O’Neel2010). The landward extent of tidal flexure is not necessarily identical to the grounding line position, as ice is thick and some of the flexure will necessarily be transferred upstream of the grounding line, but here, we treat it as a repeatable proxy for grounding line position. We allow it to move between cycles by performing a grid search around the initial estimated position where 
$\partial w/\partial x = 0$, which we first estimate as described below in the ICESat-2 grounding zone deflection data section. We repeat the optimization described above at discrete points around the initial estimated 
$x_F$ up to 1 km away in both along-track directions, at intervals of 100 m.
This implicitly allows for tidal migration of the grounding line, as well as includes any flexure of the bed that occurs below the observable surface flexure. This method is akin to a ‘free/stiff’ configuration in the parlance of Sayag and Worster (Reference Sayag and Worster2013), where the ice shelf rests on a stiff bed but the grounding line position is allowed to vary at different points in the tidal cycle. Note that we do not assume nor require the bed to be stiff, but rather that we use the combined flexure signal to infer E*.
We show the results of the boundary position grid search at each site, beam and cycle in Supplementary Material S1, and use the least-squares minimum-misfit E* and 
$x_F$ as the result for each track. We also use the shape and stability of the grid search misfit surfaces as diagnostic criteria for acceptability of the inversion (see Supplementary Material S1).
ICESat-2 flexure zone deflection data
Tidal flexure of ice shelves can be readily observed by making repeated measurements of ice shelf surface height at different points in the tidal cycle and differencing those measurements. This has been accomplished through a variety of methods, such as ground-based surveys (Vaughan, Reference Vaughan1995), satellite laser altimetry (Fricker and Padman, Reference Fricker and Padman2006; Brunt and others, Reference Brunt, Fricker, Padman, Scambos and O’Neel2010; Li and others, Reference Li, Dawson, Chuter and Bamber2022; Freer and others, Reference Freer, Marsh, Hogg, Fricker and Padman2023), radar interferometry (Rosier and others, Reference Rosier, Marsh, Rack, Gudmundsson, Wild and Ryan2017; Wild and others, Reference Wild, Marsh and Rack2018) or as here, with ICESat-2, a laser altimeter launched in 2018 with much denser spatial resolution than ICESat. We identify ICESat-2 reference ground tracks (RGTs) that cross the grounding line in three regions of the RIS: Marie Byrd Land (MBL), MacAyeal Ice Stream (MIS) and the Siple Coast (SC) (Fig. 2). We use ATL06 land ice height data (Smith and others, Reference Smith2023a) and use OpenAltimetry (J and others, Reference Khalsa2020) to identify ground tracks of interest that were approximately perpendicular to the grounding line in the NASA MEaSUREs Antarctic grounding line dataset (Mouginot and others, Reference Mouginot, Scheuchl and Rignot2017) and away from confining topography that might influence the flexure.
Repeat track measurements occur at 91 day intervals (one satellite cycle). After excluding passes from 2018, when satellite pointing was not yet optimized, as well as cycles with obvious discontinuities due to cloudy conditions or other anomalous effects, between four and eight usable passes remain for a given beam pair, which consists of a strong and a weak beam. Following ICESat-2 algorithm theoretical basis document terminology (Neumann and others, Reference Neumann2023; Smith and others, Reference Smith2023b), we use all three beam pairs gt1r/gt1l, gt2r/gt2l and gt3r/gt3l on ICESat-2, which are spaced 3.3 km apart, finding a value for E* for each cycle, for all beam pairs independently to explore the kilometer-scale variability in E*. We report the inferred E* value for each site, cycle and beam in Table 1 and the mean across all cycles at each beam at each site in Fig. 2b.
To construct the vertical tidal flexure from the repeat track ice elevation data, after processing and filtering the data as described in Supplementary Material S1, we calculate the per-cycle surface height anomaly. At each along-track point, the mean height over all N cycles is
\begin{equation}
\overline{s}(x) = \frac{1}{N} \sum_{k=1}^{N}s(x,t_k),
\end{equation}which we then subtract from the height 
$s(x,t)$ measured at each cycle to form the surface height anomaly:
\begin{equation}
w(x,t) = s(x,t) - \overline{s}(x),
\end{equation}which we interpret directly as 
$w$, the vertical tidal deflection of ice. We interpret the height anomaly as the vertical tidal deflection 
$w$ under the assumption that with enough cycles (i.e., repeat track observations), the mean height approximates the unstressed state with zero tidal amplitude. While a number of different approaches have been used in contemporaneous work to extract the tidal flexure signal from altimetry data, for instance, subtracting instead the height at the lowest measured point in the tidal cycle, or one close to the mean tidal height (Freer and others, Reference Freer, Marsh, Hogg, Fricker and Padman2023), this choice centers the spatial variability rather than the temporal variability in tidal flexure. As we make no attempt to comment here on the specific processes that might affect effective Young’s modulus at a given point in time and space, we argue that using the height anomaly is appropriate for this study.
The approximate locations of the boundaries of the flexure zone are selected visually from the ICESat-2 grounding zone deflection data: upstream of the landward limit of tidal flexure, 
$x_F$, the vertical tidal deflection should be zero, and downstream of the seaward limit of tidal flexure, it should be constant-valued. This is straightforward to identify manually for a small number of cases, and difficult to do generally (Li and others, Reference Li, Dawson, Chuter and Bamber2022). Here, because we then vary 
$x_F$ by 1 km landward and seaward of its estimated initial position, and because the exact location of the seaward boundary does not sensitively affect the result of the inversion, so long as it is far enough away to not artificially shorten the flexure zone, visual approximation of the boundary positions is sufficient.
The value of the far-field sea level, 
$A_0$, is also determined from the data as the deflection at the seaward edge of the grounding zone minus the deflection at the grounding line. This centering removes spatially uniform sea-surface anomalies, so we include no dynamic atmospheric correction. We assume 
$A_0$ is constant because ice shelf cavity bathymetry and tidal amplitude close to the grounding line are not yet well-constrained (Le Merle and others, Reference Le Merle2024). When 
$A_0$ is small, little information about ice flexure can be gleaned by the inversion and these cases are left out of our analysis.
The apparent ice surface height will change as it snows. Most products of Antarctic precipitation in this region are available only at 25 km resolution (Van Wessem and others, Reference Van Wessem2018), so this could only be modeled as a constant valued offset to the observed surface. While this is visible in some of our data as an approximately constant valued offset in height anomaly landward of the grounding line, since we effectively zero the landward boundary in our model by making 
$A_0$ the difference in deflection between the seaward and landward edge of the model, we make no further correction for snow. Snowfall on the RIS varies by site but is typically on the order of several centimeters to a meter of snow per year (Cohen and Dean, Reference Cohen and Dean2013), while the tidal amplitude is approximately 1 m about once per day (Padman and others, Reference Padman, Siegfried and Fricker2018).
Results
Across three beam pairs at three sites and ICESat-2 cycles from 2019 to 2025, we find a mean effective Young’s modulus E* = 4.7 
$\pm$ 2.4 GPa from 47 inversions. Results for each included cycle at all sites are shown in Table 1. The mean E* across all cycles at each beam pair ground track is shown in Fig. 2b, demonstrating the spatial variability in E* at kilometer scales. The modeled versus observed flexure at all three beam pairs at each site is shown in Fig. 3. Results for each track are individually plotted in Supplementary Material S1, where excluded cycles are plotted as well.
Figure 3. Observed tidal flexure at three study sites along ICESat-2 ground tracks. Each curve 
$A_{1-8}$ represents a measurement at a different point in the tidal cycle and corresponds to a unique date of observation listed in Supplementary Material S1, in its own color. Vertical flexure is measured from repeat-track surface elevation anomalies derived from ICESat-2 data. Dotted lines depict data. Solid lines depict modeled flexure. Unpaired dotted lines depict tracks that could be modeled at some but not all beam pairs.
Table 1. Inferred effective Young’s modulus (E*) in GPa at each beam and cycle used across three study sites on the RIS. Tidal cycles 
$A_{1-8}$ correspond to dates listed in Supplementary Material S1 and are not identical at different sites.
The results range from 1.1 to 9.7 GPa. The median is 4.2 GPa with quartiles 
$Q_1$ = 2.9 GPa and 
$Q_3$ = 6.5 GPa (IQR = 3.6 GPa). The standard deviation about the mean of 4.7 GPa is 2.4 GPa. Site means are similar, but variability is larger at MBL. The flexure zones are between 6 and 15 km wide. We observe that the landward boundary of the flexure zone often migrates on kilometer scales on tidal timescales. There are between 2 and 8 usable (without discontinuities or other data issues) cycles per beam pair at the sites. Figure 3 shows modeled tidal flexure as solid lines and data as dashed; dashed lines without corresponding models represent cases that were not modeled due to nonconforming data. A histogram of inferred E* results is included as Fig. 4.
Figure 4. Histogram of inferred effective Young’s modulus at all sites. The four inferred E* greater than 10 GPa, shown cross-hatched, are numerical outliers and excluded from analysis (see Supplementary Material S1) but are included here to show the distribution of results.
Discussion
We present a method for estimating effective elastic rheology near the grounding zone using tidal flexure observations from ICESat-2 and show that by varying a single effective elastic parameter in the flexure zone, we can achieve fits to flexure data that are close to observational error in some cases, with ice-penetrating radar ice thickness that varies in space from the ROSETTA-Ice campaign. Our results suggest that effective rheology near the grounding line may vary significantly in space.
Our model is linear elastic. Linear elasticity is an adequate approximation when the forcing timescale is short compared to the viscous relaxation timescale of ice in the grounding zone—that is, when the Deborah number De = 
$\frac{\tau_{ve}}{\tau_{forcing}}$ is large. Taking the viscous relaxation timescale 
$\tau_{ve}$ = 2
$\mu E$ (Turcotte and Schubert, Reference Turcotte and Schubert2014), where 
$\mu$ is the dynamic viscosity, approximated as 
$6 \times 10^{14}$ Pa s from Ranganathan and Minchew (Reference Ranganathan and Minchew2024), using E* = 4.7 GPa, the mean value from results found here, and taking the forcing (tidal) timescale on the RIS to be approximately 12 hours (Padman and others, Reference Padman, Siegfried and Fricker2018), we find De 
$\approx O(10)$. With E = 9 GPa, De is closer to 3; with E = 1 GPa it is closer to 30. This indicates that while the linear elastic approximation may be valid in some grounding zones around Antarctica, care is warranted in selection of those areas, as the effective rheology of ice is increasingly observed to have strong spatial variability due to variations in stresses, ice temperature, crystallographic fabric and grain size, interstitial liquid water content, impurities, and fractures and damage, and there may be places where the viscous timescale is comparable to the tidal timescale. The sub-basin level spatial variation in effective Young’s modulus inferred here may also have implications for flexure in tidal models that assume effective Young’s modulus is uniform (Greene and others, Reference Greene, Erofeeva, Padman, Howard, Sutterley and Egbert2024).
There is no indication that the available data support differentiation between the ice and bed elasticity. Using a single effective elastic parameter is a convenient choice that allows the uncertainty associated with doing so to be quantified. In any tidal flexure model, there is a tradeoff in uncertainty in the elastic parameters of the ice or the bed or the tide selected for use in the model, which is evident in the spatial variation we find in effective Young’s modulus. In addition to the known temperature and frequency dependence of E, which implicitly accounts for some anelastic effects, and the increasingly evident variation in the viscous properties of ice, there is no consensus yet about how the elastic properties of ice might vary as it exists in situ. There is basal and surface crevassing as well as microcracks, interstitial inclusions, larger grain sizes than laboratory experiments allow, air bubbles and other effects that might be collectively parameterized as damage.
Overall, our mean result of E* = 4.7 
$\pm$ 2.4 GPa suggests that until we can better disentangle the intersecting processes in the grounding zone, we must accept high uncertainty in effective rheology. Any attempt at parameterization of the grounding zone may require a piecewise approach with different assumptions in different locations. The magnitude of this uncertainty is similar to the result of Rosier and others (Reference Rosier, Marsh, Rack, Gudmundsson, Wild and Ryan2017), who find a factor of two of uncertainty in effective Young’s modulus. Laboratory experiments that can more closely emulate in situ ice with inclusions or damage, comparable grain sizes, and at forcing frequencies close to the tidal range would be valuable in understanding the natural laboratory of the Antarctic Ice Sheet.
We did not find obvious differences in effective Young’s modulus at high versus low tide at any site. This is in contrast to Sayag and Worster (Reference Sayag and Worster2013), who find higher elastic moduli at high tide than at low tide and argue that subglacial lubrication is enhanced at high tides. We cannot rule out due to the small number of sites that tidal effects may exist, and indeed agree that subglacial hydrology may affect observable tidal flexure. However, we argue that the elastic properties of the bed and subglacial till are effectively unknown and cannot be distinguished from the flexure of the ice, and are thus incorporated into our inference of E*. First principles modeling of the subglacial environment and how it is affected by tides, and how this translates to surface flexure, will help step toward a more accurate representation of the grounding zone in larger ice models.
A one-dimensional model minimizes the number of free parameters and makes an ill-posed problem constrained by noisy data tractable. However, if flexure gradients are steep in across-track directions, a plate bending model that incorporates flexure in both horizontal dimensions would be more appropriate. In the cases presented here, the along-track gradients are generally about an order of magnitude greater than across-track gradients, so we conclude that a 1D approximation is valid at these sites. Beam bending models, or plate bending models in 2D or 3D, with different boundary conditions, may be able to accommodate a wider subset of ice geometries than presented here, as in Ultee and others (Reference Ultee, Meyer and Minchew2020).
Out of 30 candidate RGTs around the RIS, only 3 sites were suitable enough for inclusion in the results. We select and use tracks whose flexure profiles suggest that the linear elastic approximation is a reasonable one: with smooth flexure between a flat upstream region and a flat downstream region that varies with time. We pick tracks that are roughly perpendicular to the grounding line and away from confining topography that might induce flexure in cross-track dimensions. No sites in the Transantarctic Mountains met these criteria: those tracks had different patterns of tidal deflection. Those grounding lines may also be substantially three-dimensional as the ice flows off steep topography.
Future work that seeks to parameterize grounding zone processes in larger models may need to take into account the inherent spatial variability in effective grounding zone rheology in order to accurately translate from tidal timescale processes to long-term ice sheet dynamics. Finding E* in more continuous regions, especially where across-track gradients might be significant, will help clarify the true spatial variance in effective Young’s modulus and indicate where a plate-bending approach would be more appropriate.
The choice to parameterize the grounding zone in terms of an effective elastic rheology is both motivated by the need for an independent estimation of the ice thickness and thickness gradient there in order to calculate the basal melt rate of ice near the grounding line and also justified by the uncertainty about basal conditions in the grounding zone and of the rheology of ice. The ice flexure reveals that there is evident variability in effective rheology on sub-catchment scales; what causes this variability and how it can be parameterized in large ice models will be increasingly constrained by fundamental physics models constrained by the remote sensing data record that accumulates day by day.
Conclusion
We infer an effective Young’s modulus of 4.7 
$\pm$ 2.4 GPa from ICESat-2 data by modeling the flexure zone as an elastic beam of variable thickness. Surface flexure measurements do not readily allow us to distinguish between deformation of the bed and deformation of the ice. Therefore, the inferred elastic modulus reflects the combined flexural response of both ice and bed, along with any additive tidal or subglacial hydrological effects, as well as the in situ relationship between stress and strain in ice, which is time-dependent due to its viscoelasticity. Understanding how the effective Young’s modulus varies in space at kilometer scales offers new insight into the mechanical character of the ice–bed–ocean interface.
Making models of Antarctica that accurately reflect its current state and recent past requires careful modeling of the grounding zone and better parameterizations of the processes that occur there on intersecting spatial and temporal scales. Our work shows that a linear elastic model of tidal flexure could be one method for doing so, but also that it would likely need to be part of a suite of solutions encompassing viscous, viscoelastic and viscoplastic flexure as well. The geometry of the grounding zone is complex and dynamic: different bending geometries that fit the topography on the scales over which it varies could be employed. Extending our elastic beam bending method to a viscoelastic plate bending flexure method is one clear way that more complex geometries could be represented, but the extent to which any of these potential methodologies are differentiable from one another within a small, noisy, observable signal is not yet well constrained.
The results presented here add to the constellation of estimates of the effective rheology of ice from observational and experimental methods. We are able to consider questions of variability of effective Young’s modulus on small scales only because this is one of the areas where direct ice-penetrating radar measurements of ice thickness exist. In other places, ice thickness in the grounding zone is estimated by assuming ice shelves are fully free-floating and their thickness can be calculated by measuring the height of the ice above flotation. In the grounding zone, this would bias ice thickness estimates. In order to make an independent estimation of the ice thickness near the grounding line from the surface flexure, we must first make a reasonable estimation of the effective Young’s modulus. More direct measurements of ice thickness in grounding zones would be a rich resource for the glaciological community as we attempt to untangle the intersecting physical processes in grounding zones, and are deeply pertinent to IPCC goals toward reducing uncertainty about future sea level rise.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2026.10137.
Acknowledgements
The authors wish to thank Bryan Riel for providing early assistance with data processing, Rina Cao for helping expand the scope of the work and two anonymous reviewers for their contributions. FE received funding through the NASA FINESST award 80NSSC21K1619, the John W. Jarve (1978) Seed Fund at MIT and the NEC Corporation Fund for Research in Computers and Communications at MIT. BM received funding through NSF-NERC award 1853918, the John W. Jarve (1978) Seed Fund at MIT and the NEC Corporation Fund for Research in Computers and Communications at MIT. CRM was supported by the Army Research Office (78811EG).
Competing interests
BM and CM are co-founders of Arête Glacier Initiative (areteglaciers.org), where we maintain affiliations through our allowance for outside professional activities provided by our respective universities. Arête is a non-profit organization (currently a fiscally sponsored organization of the 401(c) 3 Digital Harbor Foundation) founded in 2024 to provide funding to glaciological research focused on sea-level rise. BM, FE and CM declare no financial, personal or professional conflicts of interest.
Open access
