Ocean Coupling Limits Rupture Velocity of Fastest Observed Ice Shelf Rift Propagation Event

The Antarctic ice sheet is buttressed by floating ice shelves that calve icebergs along large fractures called rifts. Despite the significant influence exerted by rifting on ice shelf geometry and buttressing, the scarcity of in situ observations of rift propagation contributes considerable uncertainty to understanding rift dynamics. Here, we report the first‐ever seismic recording of a multiple‐kilometer rift propagation event. Remote sensing and seismic recordings reveal that a rift in the Pine Island Glacier Ice Shelf extended 10.53 km at a speed of 35.1 m/s, the fastest known ice fracture at this scale. We simulate ocean‐coupled rift propagation and find that the dynamics of water flow within the rift limit the propagation rate, resulting in rupture two orders of magnitude slower than typically predicted for brittle fracture. Using seismic recordings of the elastic waves generated during rift propagation, we estimate that ocean water flows into the rift at a rate of at least 2,300 m3/s during rift propagation and causes mixing in the subshelf cavity. Our observations support the hypotheses that large ice shelf rift propagation events are brittle, hydrodynamically limited, and exhibit sensitive coupling with the surrounding ocean.

Rifting remains a challenging physical process to observe.Many rifts in Antarctica's largest ice shelves, like Ross Ice Shelf, appear to be stable on decadal timescales with minimal propagation in the observational record (Walker et al., 2013).In contrast, rifts on highly dynamic ice shelves like Pine Island Glacier initiate, propagate, and calve icebergs every few years (Jeong et al., 2016;Olinger et al., 2019;Walker et al., 2013).While remote sensing provides an accurate and effective method of measuring many aspects of ice shelf evolution, the wide range of rift propagation timescales prevents the full spectrum of rift behavior from being observed by remote sensing alone.Because the revisit period of most satellite instruments is several days apart, remotely-sensed observations cannot distinguish between rift propagation that occurs continuously over a timescale of several days and rift propagation that occurs rapidly over a timescale of hours or seconds.The need for high-resolution observations of rift propagation has been answered in part by deploying seismic arrays to continuously monitor the elastic wave emissions from rifts on short timescales (Bassis et al., 2005;Olinger et al., 2019).However, despite the promising results of rift seismology, the logistical complexity and hazard of field campaigns in active areas of ice shelf deformation mean that only a handful of seismic arrays have been deployed near rifts.Furthermore, seismic studies have been unable to capture any instances of truly rapid rifting despite evidence that such events do occur (Banwell et al., 2017).
In this study, we analyze rift propagation at Pine Island Glacier Ice Shelf (PIG), a fast-flowing ice shelf in West Antarctica that was the single largest Antarctic contributor to sea level rise in the period 1979-2017 (Rignot et al., 2019).Since 1992, tabular icebergs have calved from PIG every 2-6 years along rifts that propagated from the northern and southern shear margins, maintaining a relatively consistent ice front position and orientation (Arndt et al., 2018).In 2015, calving occurred along a rift that initiated in the ice shelf's center for the first time, resulting in substantial ice front retreat and reorientation (Jeong et al., 2016) that has continued to the present.Before this change in ice front geometry, the last calving event occurred in 2013 along a rift that propagated from the northern shear margin across the ice shelf, hereinafter referred to as R2011 for the year of its initiation.Here, we overcome the perennial limitation of studies of calving processes -namely, a lack of in situ observations (Benn et al., 2007)-and present the first-known near-field seismic observations of a large ice shelf rift propagation event.

Identifying Rift Propagation in SAR Data
We manually examine synthetic aperture radar imagery collected by the TerraSAR-X (TSX) satellite (Pitz & Miller, 2010) to identify an episode of rapid propagation during the rifting that preceded the 2013 calving event.One of the northern shear margin fractures, hereinafter referred to as R2012, propagated across the ice shelf and connected with R2011.TSX imagery from 8 May 2012 04:04 UTC shows the PIG ice shelf before the episode of rift extension (Figure 1a).The rift R2011 spanned 33.79 km across the ice shelf from the northern shear margin, and a band of ∼20 shorter parallel fractures spanned ∼5 km from the northern shear margin.The data show no major rifts besides R2011.TSX imagery from 11 May 2012 03:13 shows the PIG ice shelf after the episode of rift extension (Figure 1b), providing a 3-day time window around the episode of rift extension.A high-resolution digital elevation model of PIG (Shean et al., 2019) shows that the tip of R2012 was located in a basal trough before propagation, consistent with previous observations that suggest basal channels strongly influence rift propagation on PIG and other ice shelves (Alley et al., 2019;Dow et al., 2018).

Identifying Rift Propagation in Seismic Data
We identify the rift seismic signal within data recorded by three Nanometrics Trillium 120 seismometers deployed on the PIG ice shelf and three Nanometrics Trillium 240 seismometers deployed across West Antarctica (D. Holland & Bindschadler, 2012).After correcting for the instrumental response, seismic data recorded in the time window established by TSX imagery contains a single notable signal, recorded on 9 May 2012 at 18:03 (Figure 1c).We hypothesize that this signal, the largest amplitude signal within the three-day window, was generated by the extension of R2012 observed in TSX imagery.We use the cumulative amplitude distribution of the signal to estimate a duration of 2.09 hr, longer in duration than all other signals in the time window by an order of magnitude.The May 9 event has a peak vertical ground velocity of 0.234 mm/s and peak vertical ground displacement of 0.195 mm at a distance of 12 km and uniquely contains significant energy at periods up to 1,000 s.Despite correcting for the instrumental response, the sensitivity of the Trillium 120 seismometer is reduced below its natural period of 120 s, suggesting that the amplitudes recorded between periods of 120 and 1,000 s may underestimate actual ice shelf velocities at those periods.Between 1,000 s and 1 s periods, the signal exhibits high-frequency-first dispersion characteristic of flexural gravity (FG) waves (Figure 1d), a wave type that propagates as a coupled beam flexural and ocean surface wave (Abrahams et al., 2022;Press & Ewing, 1951;Sergienko, 2017;Squire, 2007).At frequencies above 1 Hz, the signal consists of P-waves and Rayleigh waves that gradually increase in amplitude before abruptly decaying after 300 s (Figure 1e).These higher-frequency phases are also recorded by regional POLENET stations DNTW, THUR, and UPTW, respectively located 250, 294, and 360 km from PIG.

Mapping Rift Extent
We employ a semi-automated scheme to identify the extent of R2012 before and after propagation.We use TSX images from 5 May 2012 03:22:11 and 11 May 2012 03:13:39, which were captured from similar incidence angles and span the same spatial extent.To remove the effect of ice shelf advection, we cross-correlate windows containing the rift tip from each TSX image to obtain the optimal shift between the two images.We then use the computed shift to align the two images.To increase the efficiency of image processing, we then downsample both images by a factor of 10.To measure the increase in length of R2012, we normalize the images from May 5 and May 11 such that pixels with values close to 1 correspond to dark features like rifts.We then subtract the pre-extension image from the post-extension image to remove all features constant between May 5 and May 11, including shear margin fractures and R2011.We extract the largest 1-valued region from the differenced images, corresponding to the increase in the area of R2012.We then skeletonize the binary rift image to obtain a thinned representation of the rift's path with a width of one pixel (Maragos & Schafer, 1986).We then measure the length of the skeleton's main branch in pixels and multiply by the TSX image's pixel size to extract the increase in length between May 5 and May 11.Finally, we sum the binary rift image to obtain the area of the rift in pixels, multiply by the TSX image's pixel size to obtain the rift area in square meters, then divide by the increase in length of R2012 to obtain an estimate of the average rift width.We estimate an increase in length of 10,530 m and a final average width of 150 m.We follow the same procedure to estimate the initial length and width of R2012, finding an initial length of 3,890 m and an initial width of 90 m.Based on the TSX images' pixel size of 1.66 m and the downsampling performed in our workflow, quantities estimated using TSX images have an uncertainty of ±16.6 m.

Seismic Location
To locate the candidate rift event, we first employ a grid search algorithm using arrival times at locally and regionally-deployed stations.To obtain the relative arrival times of high-frequency waves (1.5-5 Hz) at each station, we cross-correlate the filtered signal recorded at PIG2, the closest station to R2012, with the filtered signal recorded at each other station.Propagation paths from PIG to regional stations contain transitions between floating ice and grounded ice, with significant variation in ice thickness and topography.We opt to use estimated phase velocities, which represent the phase velocity averaged over the propagation path, instead of a single phase velocity from literature due to the lateral heterogeneity in velocity structure.We estimate the phase velocity of the waves recorded at each regional station by dividing the known distance between PIG2 and each station by the difference between the arrival time at PIG2 and each station.We carry out this procedure using vertical, north-south, and east-west component data recorded at regional stations THUR, BEAR, DNTW, and UPTW.For local stations PIG2, PIG4, and PIG5, we assume a constant phase velocity of 2,000 m/s corresponding to the Rayleigh wave speed in ice (Kim et al., 2010).We then conduct the grid search by iterating through possible origin times and spatial locations and computing the expected arrival time at each station using the corresponding phase velocity.We calculate the root mean square error (RMSE) between the observed and expected arrival times for all components and stations, giving a single estimate of the misfit in arrival times across the array.We then calculate RMSE for each possible origin time and for every spatial point in a regular grid to obtain a map of error.The event location is finally determined by identifying the spatial point and origin time that correspond to the lowest RMSE.
To further constrain the source location, we use the polarization direction of horizontal waves recorded at on-ice stations PIG2, PIG4, and PIG5 to compute an epicentral back azimuth.The recorded high frequency waves 10.1029/2023AV001023 5 of 12 (1.5-5Hz) consist of waves with back-and-forth particle motion and elliptical particle motion with a common polarization direction.Elliptical particle motion is diagnostic of Rayleigh waves.Because the phases with back-and-forth particle motion share a polarization direction with the recorded Rayleigh waves, we infer that the seismogram contains P-waves and Rayleigh waves, both wave types that are polarized in the direction of propagation.Determination of the polarization direction of recorded high frequency waves therefore reveals the direction from which they propagated to the array.By performing principal component analysis (PCA) on the east-west and north-south seismograms, we obtain the PCA first component, a vector corresponding to the direction along which the majority of the variation in the data occurs.We infer the polarization direction from the PCA first component, which corresponds to one of two possible propagation directions separated by 180°.To resolve this 180-degree ambiguity, we identify the two stations farthest from the array centroid in both possible directions of propagation, which are expected to record the first arrivals for incoming plane waves from either propagation direction.We then adjust the sign of the PCA first component to match the propagation direction whose predicted first arrival agrees with the observed first arrival.We repeat this procedure using data recorded at each station and sum the PCA first component vectors from each station to obtain an average propagation direction.Finally, we retrieve a back azimuth by taking the arctangent of the quotient of the two elements of the PCA component vector.We repeat the entire procedure for each 50 s time window in the event, resulting in a distribution of back-azimuths calculated for each time window within the event.We obtain a single event back azimuth by taking the circular mean of the back-azimuths calculated from each time window, with the back azimuth from each time window weighted by the norm of the summed PCA components across the array for that window.

Ocean-Coupled Fracture Modeling
We model the coupled ocean-rift system using simple linear elastic fracture mechanics and fluid dynamics.Given the relatively limited observations of rift growth, we have pursued a simplified model.Modeling efforts with increased complexity inherently lead to a greater number of model variables that, in our case, cannot be compared to observations.In the absence of more detailed observations, added complexity offers little benefit compared to the simplified case.
We utilize a coordinate system where the ice shelf base has a vertical position of 0, the ice front has a position of 0 in the y-direction, and the back of the rift has a position of 0 in the x-direction (Figure 2a).Our modeling employs the following equations, which are derived or discussed in more detail in Text S1-S6 of Supporting Information S1.
The conservation of water mass flow rate into the rift (derived in Text S1 of Supporting Information S1) is: where η is the water level within the rift, w is the width of the rift, and u is the water flow rate, H c is the subshelf cavity height (assumed uniform), and L is the rift length.Although spatial variations in H c exist (Muto et al., 2016;Shean et al., 2019) and could play a role in higher order dynamics, we assume uniformity in order to describe only the most basic features of ocean-rift coupling.
We write the conservation of fluid momentum (derived in Text S2 and S3 of Supporting Information S1) as: where L c is the horizontal position of the rift and H w is the height from the ice shelf base to the hydrostatic water line.This equation arises from a particular case of the unsteady Bernoulli flow approximation.
The depth integrated perpendicular rift opening stress (discussed in Text S4 of Supporting Information S1) is: where R xx is the ice shelf extensional stress, H i is the ice shelf thickness, ρ w is the density of water, and ρ i is the density of ice.
The rift tip equation of motion (discussed in Text S6 of Supporting Information S1) is: where c r is the Rayleigh wave speed, K I is the stress intensity factor experienced by the rift and K c is the fracture toughness of ice.
We seek to solve the system of ordinary differential equations (ODEs) defined by Equations 1, 2, and 5 for L and η.We utilize a Python implementation of the LSODA ODE solver (Hindmarsh, 1983;Petzold, 1983) that handles systems of equations with the form, that is, systems that only have a dependence on the first derivatives x′ of the state vector x.Through algebraic manipulation, we write Equations 1, 2, and 5 in this form.This requires introducing a variable χ = η′ so the second derivative η″ can be written as a first derivative.This also requires us to obtain an equation for χ″.We accomplish this by rearranging Equation 2, analytically computing all the necessary derivative terms, and substituting.This process yields a set of equations written in a convenient way for numerical solutions.However, the equations shown above are more physically interpretable, so the forms used in obtaining a numerical solution are not reproduced here.
In our simulations, we first compute the ice shelf extensional stress R xx necessary for a fracture of the prescribed initial geometry to experience a stress intensity factor equal to the fracture toughness of ice.We then add a small stress perturbation to the variable R xx compute the stress and stress intensity factors corresponding to the perturbed value of R xx assuming an initial steady-state water surface height of η = H i ρ w /ρ i .This initiates fracture propagation at a rate determined by Equation 5. R xx is held fixed through the simulation, and the stress applied to the rift evolves through time as η changes.
In reality, we have little independent constraint on the short-timescale variability of ice shelf stresses leading up to the observed propagation of R2012.Ocean swell, tides, and stick-slip motion have all been shown to create stress perturbations in the range of ∼1 kPa (Lipovsky, 2018;Olinger et al., 2019;Wiens et al., 2022), and a combination of these effects may have provided the increase in rift stresses that initiated the propagation of R2012.The sensitivity of modeled rupture velocities to ice shelf geometry and R xx are shown in Figures S2-S4 of Supporting Information S1.We find that the rupture velocity predicted by our model is highly sensitive to the exact magnitude of ice shelf extensional stress R xx that is applied to initiate rift propagation.This sensitivity suggests that observed rupture velocities during rifting may be used to understand the short term stresses in ice shelves, in good analogy with the analysis of tectonic earthquake ruptures (Dunham, 2007).

Origin of Observed Seismic Signal
A previous study of fracture at PIG identified FG waves generated by gradual propagation of R2011 (Olinger et al., 2022), suggesting the extension of R2012 is a reasonable source for the May 9 event.In addition to rift propagation, FG waves on ice shelves are generated by incoming ocean waves (Chen et al., 2020).However, ocean wave sources cannot account for seismic phases recorded by regional stations hundreds of kilometers away, and we therefore conclude that incoming ocean waves must not have generated the rift signal.Additionally, the spectrum of the May 9 event is markedly different than teleseismic earthquake spectra recorded by the same instrument (Figure S1 in Supporting Information S1), and we conclude that teleseismic waves must not have generated the May 9 event signal.A grid-search inversion of arrival times at locally and regionally-deployed seismic stations finds its lowest-error region where R2012 connects to R2011 (Figure 1b), further supporting the hypothesis that the extension of R2012 generated the May 9 event.The polarization of waves recorded at locally-deployed stations corresponds to a back azimuth of 308.1 ± 6.2° (Figure 1b), in agreement with the back azimuth of R2012, confirming that the recorded waves propagated to the local seismic stations from the direction of R2012.Because the best-fit event location coincides with R2012 and because both teleseismic and ocean wave sources are inconsistent with the seismic observations, we conclude that the May 9, 18:03 event was the seismic signal generated by the extension of R2012.

Rate of Observed Rift Propagation
To understand the dynamics of the observed rift propagation, we estimate the rupture velocity using the duration of radiated body and surface waves and the increase in length estimated from TSX imagery.Radiated body and surface waves gradually crescendo, consistent with an accelerating rupture or a seismic source that moves progressively closer to the seismometers (our limited station geometry precludes distinguishing between these two scenarios).Seismic waves then abruptly stop after 300 s (Figure 1e), indicating the conclusion of propagation when R2012 collided with R2011.Such a "stopping phase" is highly unusual; stopping phases are not typically observed in tectonic earthquakes, for example, We thus infer that the observed 10.53 km of rift extension occurred over 300 s, corresponding to an average rupture velocity of 35.1 m/s.To the best of our knowledge, this is the fastest rift propagation speed ever observed.R2012 extended over a duration of time two orders of magnitude below the Maxwell time of ice, which is around 11 hr (Ultee et al., 2020), supporting the hypothesis that the observed rift extension occurred through dynamic brittle fracturing.However, elastodynamic theory predicts opening mode fracture propagation at rates approaching the Rayleigh wave speed of the fracturing material (Freund, 1990), which is between 1,500 and 2,000 m/s in the ice (Kim et al., 2010).Why, then, did R2012 propagate two orders of magnitude below the Rayleigh wave speed?

Modeled Ocean-Coupled Rift Dynamics
We hypothesize that coupling between rift propagation and water flow within the rift explains rupture at a small fraction of the Rayleigh wave speed.To test the hypothesis, we develop a simple model of rift propagation that couples brittle fracture and water flow.The rift is represented by a sharp fracture in a linear elastic plate subject to uniform far-field tension.Water flows through the subshelf cavity into the rift is represented using the unsteady Bernoulli flow approximation.To initiate rift extension, we apply the stress required for a fracture of a chosen initial geometry to experience a stress intensity factor that just exceeds the fracture toughness of ice.As the rift propagates, the total volume within the rift increases, and water rushes in to fill the rift.However, water flow into the rift is not rapid enough to maintain the hydrostatic water line, causing a reduction in the average water height within the rift and a decrease in the depth-integrated water pressure acting to open the rift (Figure 2a).The lower water pressure reduces the total resolved stress that drives rift opening, and in turn, limits the rate of propagation.This effect results in a far lower rupture velocity than predicted for a rift with a static water height.If propagation stops abruptly, fluid inflow continues due to inertia and overshoots the steady-state water line, resulting in simple harmonic oscillations about the steady-state water line (e.g., Figure 2b after t = 300 s).These simple harmonic oscillations occur at the sloshing period,   = 2 √ ∕ , where M = L c w/H c + H w − H c is a measure of the effective cross-sectional area of water being transported, L c is the distance from the rift to the ice front, w is the width of the rift, H c is the distance from the seafloor to the base of the ice shelf, and H w is the distance from the seafloor to the water surface.
To test whether the proposed model can explain the observed propagation rate of R2012, we model R2012 using an initial length of 3.89 km and an initial width of 90 m that were measured from TSX imagery.We assume that the ice shelf has a uniform ice thickness of 400 m, estimated from a high-resolution digital elevation model of PIG from 2012 (Shean et al., 2019), and a uniform water depth of 840 m, estimated from a gravity-derived model of Pine Island Bay bathymetry (Muto et al., 2016).We then obtain the magnitude of R xx by computing the stress required for a rift with the measured initial geometry to begin unstable propagation and identifying the additional stress needed for the rift to grow 10.53 km in 300 s.The modeled rift begins propagating when an extensional stress of approximately 161 kPa is applied.This is consistent with a previous estimate which found that the central region of PIG ice shelf had a mean R xx of approximately 124 kPa (Lai et al., 2020) (Figure S5 in Supporting Information S1).Once propagation begins, the rift propagates 10.53 km over a duration of 300 s at an average rate of 35.1 m/s (Figure 2e), in agreement with the observed rupture velocity.The ability of the model to reproduce the observed rupture velocity strongly supports the hypothesis that coupling between rift propagation and water flow limited the observed rupture velocity to a small fraction of the Rayleigh wave speed.Without fluid coupling, a fracture of the same initial geometry subjected to the same magnitude of stress reaches 99% of the Rayleigh wave speed in 300 s, highlighting how significantly the mechanism we propose influences the dynamics of rift propagation.

Discussion
Our observations suggest that ice shelf rift propagation occurs more rapidly than suggested by previous seismic and remotely sensed observations (Bassis et al., 2005;Olinger et al., 2022;Walker et al., 2013).An immediate conclusion from this observation is that the timescale of fracture at R2012 is well within the regime of brittle fracture.Yet our observations and modeling show that the rift propagation of R2012 was not purely governed by the laws of linear elastic fracture mechanics because rift propagation was slowed down by interaction with the ocean.Several aspects of R2012 warrant further discussion.

Hydrodynamically-Controlled Fracture
Water has long been thought to favor fracture growth in glaciers.The presence of water enables vertical crevasses to propagate through the entire ice thickness of glaciers and ice shelves (van der Veen, 1998; Weertman, 1973) that may result in calving depending on water supply (Zarrinderakht et al., 2022).Meltwater accumulation in supraglacial lakes can load ice shelves sufficiently to generate flexural stresses that trigger widespread ice shelf fracture and collapse (Banwell et al., 2013), wherein increased meltwater supply encourages more rapid and pervasive fracturing (Robel & Banwell, 2019).Water pressure also makes a significant contribution to the stresses that drive rift propagation, and rifts filled with mélange, a semi-consolidated mix of ice and snow, appear more stable than rifts filled with open ocean water (Huth et al., 2023;Larour et al., 2021).In each case, water has a primary role as a destabilizing mechanism for fracture.Our observations and models, however, show that water has a dual role in rifting.While the stress exerted by water within rifts is necessary to initiate propagation, the dynamics of water flow during propagation serve to limit the rupture velocity of the rift.Whereas the largest fracturing events on land, that is, tectonic earthquakes, are ultimately inertially limited (Dunham, 2007), our observations and models imply that the largest fracturing events in ice, for example, ice shelf rift propagation events, are ultimately hydrodynamically limited.Our results suggest that, especially for cases of rapid fracture on short timescales, attempts to model the propagation of fluid-filled fractures should not neglect the dynamics of water flow.While the implications of our model most readily apply to ice fracture, we note that the dynamic role of fluid flow in controlling rift propagation shows some similarity to the role of fluid in other geophysical settings such as growth of pressurized volcanic dykes (McLeod & Tait, 1999) and earthquake slip along pressurized faults (Zhu et al., 2020).

Rapid Rifting and Ice Shelf Stability
Despite the slowing effect of rift hydrodynamics, the propagation of R2012 is nonetheless the fastest observed episode of rift propagation.Rapid rift propagation represents a possible mechanism of sudden ice front retreat or ice shelf collapse.While previous observations of rift propagation have overwhelmingly captured gradual or episodic propagation (Bassis et al., 2005(Bassis et al., , 2007;;Jeong et al., 2016;Walker et al., 2013), we show that rift propagation on the order of 10 km can occur in a matter of minutes.It is unknown whether this represents a rare class of rift behavior or a relatively common class of rift behavior that has remained undetected until now due to the low temporal resolution of remotely-sensed observations and a scarcity of ice shelf seismic deployments.As PIG continues to accelerate, elevated stresses and shear margin weakening are expected to enhance rift propagation (Lhermitte et al., 2020;Lipovsky, 2020), which can, in turn, lead to ice front retreat, buttressing loss and further acceleration (Joughin et al., 2021).Our observations suggest that such feedback at PIG and other unstable ice shelves across Antarctica may progress more rapidly than anticipated.

FG Waves Generated by Rift Propagation
We use the wave impedance tensor (Lipovsky, 2018) to compute the maximum flexural stresses carried by FG waves recorded at each local station and estimate a mean flexural stress of 3.26 kPa, consistent with typical ocean wave-induced flexural stresses on Ross Ice Shelf (Aster et al., 2021;Lipovsky, 2018) and potentially large enough to trigger additional fracturing within the ice shelf.However, PIG typically experiences a lower degree of ocean wave excitation than Ross Ice Shelf (Chen et al., 2018), so FG waves generated by rift propagation may exert a greater influence than ocean waves on the stability of fractures on PIG.
The 600 s dominant period of the recorded FG waves is between the greatest ice shelf resonance period (∼1,600 s) and the sloshing period (T slosh ∼ 100 s), suggesting that both of these processes are involved in generating the observed FG wavefield.Accounting for radiative losses at the ice front using the relevant reflection coefficient (Abrahams et al., 2022) results in an underestimate of the e-folding duration (i.e., time to achieve decay by a factor of 1/e) of the wavefield as 16.7 min.However, it takes 38.2 min for recorded FG waves to decay by a factor of 1/e.The e-folding duration is therefore plausibly attributed to wave generation by water sloshing within the rift that continues after rift propagation ceases (see, e.g., Figure 2b).

Mixing Induced by Rapid Rifting
We infer that large rift propagation events induce diapycnal mixing in the subshelf cavity, that is, vertical flow in the presence of horizontal density surfaces (P.R. Holland et al., 2019;Jacobs et al., 2011).In the context of smaller-scale calving, Meredith et al. (2022) recently observed such mixing following a calving event with potential energy change (0.6-2.4) × 10 12 J.At R2012, we estimate the potential energy change 1.2 × 10 12 J over the 5-min duration of rift propagation and 1 × 10 15 J over the subsequent several days of rift opening (Text S7 of Supporting Information S1).Whether the internal tsunami mixing mechanism proposed by Meredith et al. (2022) is able to operate at the timescale of the longer-duration potential energy change depends on the water column stratification through the buoyancy oscillation frequency (Gill, 1982), information which is not available during the time of R2012.However, for both cases, the scale of energy associated with vertical diapycnal flow implies significant subshelf mixing during rift growth, contrary to earlier reports (Meredith et al., 2022).
During the several minutes of rift propagation, we calculate the vertical water volume flux to be at least 2,300 m 3 /s.PIG is known to have a complex basal topography with pervasive longitudinal basal crevasses that penetrate as much as 30% of the ice thickness (Vaughan et al., 2012) and that are perpendicular to the propagation direction of the rift.On other ice shelves, such features have been shown to guide the direction of rift propagation (De Rydt et al., 2018).Such basal crevasses do not play a significant role in our simplified model since we do not attempt to model the rift propagation path.If the R2012 rift did follow a basal crevasse with a height of 30% the ice thickness, this would reduce our vertical flow estimate by the same percentage.
Rifting-induced mixing suggests the existence of positive feedback between these processes.Despite significant thermocline variability, sub-shelf waters in the vicinity of the rift R2012 are deep enough to consistently reach the depth of warm circumpolar deep waters (Christianson et al., 2016).Rift propagation in this setting may therefore elevate isothermal contours and cause warming of the ice-ocean interface.Localized and repeated rift propagation in areas like the northern shear margin of PIG (Figure 1a) may then initiate feedback wherein rift propagation induces mixing and localized melting that contributes to marginal weakening (Lipovsky, 2020) and the formation of basal melt channels (Alley et al., 2019;Dow et al., 2018), thereby promoting further rift propagation.

Conclusions
We conclude that rifts can propagate rapidly through brittle fracture and that the ocean exerts a profound influence on rift propagation.We observe rift propagation at a rate of 35.1 m/s, suggesting that rifting may act as a mechanism of ice shelf destabilization over short timescales.Seismic recordings show that rapid rift propagation generates FG waves that transmit flexural stresses on the order of 1 kPa across the ice shelf and may be capable of triggering subsequent fracture.We find that the dynamics of water flow during rifting can limit the rate of propagation to a small fraction of the expected rate of brittle fracture, and that rapid rift propagation induces significant water flow into the rift, which we infer drives mixing in the subshelf cavity.We therefore add to the body of literature documenting diverse ice-ocean interactions and demonstrate that extreme ice shelf sensitivity to ocean conditions extends to the fine-scale dynamics of rift propagation.

Figure 1 .
Figure 1.TSX and seismic data surrounding 10.53 km of rift propagation.All seismic data are vertical velocity seismograms recorded by station PIG2.(a) TSX image from 8 May 2012 04:04 (t 1 ) before propagation of R2012.Black triangles show on-ice seismic stations.The inset shows the location of PIG within Antarctica.(b) TSX image from 11 May 2012 03:13 (t 2 ) after propagation of R2012.Contours show the logarithm of root mean squared error (RMSE) in arrival time residuals from a grid search of possible locations of the May 9 signal computed using on-ice seismic stations and regionally-deployed seismic stations (not shown).Arrow denoted by θ shows backazimuth computed using the polarization of seismic waves recorded by on-ice seismic stations.(c) Seismogram spanning the time window between TSX images.Time period shown in D is highlighted in pink.(d) Seismogram of rift event filtered between 0.001 and 1 Hz.In this frequency band, the signal is dominated by flexural gravity waves.Resonance of ice shelf modes results in an event duration on the order of hours.Time period shown in E is highlighted in yellow.(e) Seismogram of rift event filtered above 1 Hz.In this frequency band, the signal is dominated by P waves and surface waves.The abrupt decay of the rift event signal 300 s after the onset of the event indicates the conclusion of rift propagation.

Figure 2 .
Figure 2. Ocean-coupled model of rift propagation.In all panels, black curves show a rift coupled to hydrodynamics, and orange curves show a rift uncoupled to hydrodynamics.(a) Illustration of the proposed mechanism for rift propagation at a small fraction of the Rayleigh wave speed (c r ).(b) Modeled perturbation from hydrostatic water level during rift propagation.In the coupled case, the water level initially decreases because water flow into the rift is not fast enough to maintain the hydrostatic water level.Once propagation concludes, flow overshoots the hydrostatic water level and oscillates at the sloshing period T slosh .In the uncoupled case, the water level in the rift does not change.(c) Modeled rift propagation rate through time.In the coupled case, decreasing water pressure limits propagation to a small fraction of c r .In the uncoupled case, propagation approaches c r rapidly.(d) Modeled rift length through time for 500 s of rift propagation.(e) Modeled rift length through time for 300 s of rift propagation.In the coupled case, the rift length increases by 10.53 km at an average rate of 35.1 m/s, in agreement with the seismically derived propagation rate of R2012.In the uncoupled case, the rift propagates at a rate far exceeding the observation.