Crevasse Propagation on Brittle Ice: Application to Cycloids on Europa

Existing lineaments on the surface of the Jovian moon Europa are thought to be the result of ongoing brittle crack formation in the elastic regime. Arcuate features are called cycloids and can be modeled using linear elastic fracture mechanics. Here, we build on existing terrestrial models of rift propagation and extend them to cycloids on the moon. We propose that these cracks tend to grow as a series of nearly instantaneous events, spaced by periods of inactivity. The behavior is similar to what is observed on Antarctic ice shelves, where rifts can remain dormant for years. We argue that dormant periods between growth events could explain the presence of cycloids on Europa even without invoking secular motion of the crust. Furthermore, being able to model propagation events and their timing should help future missions exploring the moon.


Introduction
Observations coming from past planetary missions, in particular the Galileo spacecraft, revealed lineaments and fractures across Europa's surface. Various fracture models and formation scenarios have been proposed, although available data are not sufficient to fully understand how the surface of the moon deforms (e.g., Rhoden et al., 2015). However, one of the main requirements for the formation of these type of cracks is the existence of liquid water underneath the icy crust (Carr et al., 1998;Pappalardo et al., 1999). Radio doppler tracking and magnetic measurements performed by Galileo revealed a water ocean below the global shell (Anderson et al., 1998;Khurana et al., 1998;Kivelson et al., 2000). Although stratification scales are still poorly constrained, the ocean should be around 100 km deep (e.g., Sohl et al., 2002), while the thickness of the frozen ice shell ranges from a few to more than 30 km (e.g., Billings & Kattenhorn, 2005). Tidal stress produced by gravitational interaction with Jupiter is the main force driving cold ice on Europa to fracture and form these crevasses (e.g., Greenberg et al., 1998;Hoppa et al., 1999aHoppa et al., , 1999b.
Some of the lineaments are characterized by arcuate and cuspate segments of around 100 km in length. These features are called cycloids, and their orientation can be explained by formation due to tensile cracking (Hoppa et al., 1999b(Hoppa et al., , 2001). In the model of Hoppa et al. (1999b), cycloidal arcs are formed in response to the constantly changing tidal stress magnitude and orientation along Europa's orbit. The stress eventually reaches the ice strength and allows the crack to evolve and follow curved patterns. The rate of change of the tidal stress can be used to derive an apparent growth rate of cycloids of roughly 3 km/hr. By constraining the growth rate, Hoppa et al. (1999b) assumed that a cycloidal arc is formed during a single orbit of the moon around Jupiter. However, low temperatures at the surface of Europa (about 100 K, Spencer et al., 1999) and high strain rates (due to the constantly rotating tidal stress) also imply that ice deforms elastically (Rist & Murrel, 1994). Therefore, cracks on brittle ice should propagate at high fractions of the speed of sound (Lee et al., 2003). The difference in magnitude between the two speeds can be explained by the fact that cracks on ice are generally evolving at nearly instantaneous rates through several events, similarly to what happens on terrestrial ice shelves. Indeed, on Earth, ice shelf rifts are formed through episodic cracking events that follows periods of inactivity, when the crack is dormant (Banwell et  • Supporting Information S1 • Figure S1 • Figure S2 • Figure S3 • Figure S4 Correspondence to: M. Poinelli, mattia.poinelli@gmail.com 2013). As consequence, the apparent growth rate, for example, measured by remote sensing, is substantially slower than the actual propagation rate. Marshall and Kattenhorn (2005) improved on the cycloid growth model of Hoppa et al. (1999b) by describing the cusp formation in response to resolved shear stress at the cusp itself, which results from tailcracking processes. Hurford et al. (2007) further updated the cycloid model with the introduction of secular stress contribution and the variation of material parameters along the cycloid. This allowed improved fits of existing cycloids by introducing a shift in longitude due to secular rotation of the shell (e.g., due to nonsynchronous rotation). Furthermore, Rhoden et al. (2010) implemented a parameter-searching algorithm that is capable of quantitatively inferring rotational parameters of Europa by matching cycloids' shape with observations of the Galileo mission. The latter model suggests that diurnal stress due to nonzero eccentricity, obliquity, and physical libration represents the tidal model that better reproduces the shape of observed cycloids. In summary, cycloids have been extensively studied although their propagation rate has not been sufficiently described in past models. We contribute to this question by introducing a fracture mechanics approach that is capable to model episodic and instantaneous cracking events.
On Earth, large availability of in situ data and observations from remote sensing can be used as a reference to better understand how ice sheets and crevasses evolve. Among numerical models of ice fractures, one of the most straightforward approaches is linear elastic fracture mechanics (LEFM; Tada et al., 2000). This type of model approximates ice as an elastic medium, allowing for discontinuities in the structure. For example, Van der Veen (1998a, 1998b adopts LEFM to calculate equilibrium lengths of surface and bottom crevasses on terrestrial glaciers. Other models that compute horizontal propagation of ice rifts in Antarctica used techniques that can be derived from LEFM theory (Hulbe et al., 2010;Larour et al., 2004b). In addition, LEFM models of crevasse propagation have been tested and validated against a series of observational campaigns in Iceland (Mottram & Benn, 2009). The LEFM approach has been previously applied to Europa with focus on the potential of cracks to penetrate through the entire ice crust, offering a direct connection between the surface and deep ocean (Crawford & Stevenson, 1988;Lee et al., 2005;Qin et al., 2007;Rudolf & Manga, 2009). The LEFM methodology could also provide insights into directions and angles of tailcrack formation at the surface (Marshall & Kattenhorn, 2005). However, no previous study has dealt with horizontal propagation of cycloids on Europa with terrestrial-based models of ice rifts. Although knowledge of surface deformation on Europa is limited, fracture mechanics tools and terrestrial analogs provide first-order estimations, which could be validated by future missions.
In the current work, we follow the LEFM approach of Larour et al. (2004b) to compute crevasse propagation rates of cracks on brittle ice. As starting point, we fix orbital and rheological parameters, and we use elastic fracture mechanics to model the episodic propagation events. Differently to past studies, which assumed that a cycloidal arc is formed continuously, our model is capable to calculate the growth rate of cycloids by including dormant periods between propagation events. Periods of inactivity between single events depend on the local stress-strain rate state and the length of the existing crack. After presenting the LEFM methodology and the results of our model, we discuss the implications that these new findings bring to our knowledge of Europa and to its future exploration.

Methods
Surface stresses on Europa can be calculated through the computation of the tidal potential for internally differentiated spherical bodies (Harada & Kurita, 2006;Jara-Orué & Vermeersen, 2011;Wahr et al., 2009). Tidal stresses on the ice layer change in magnitude and orientation as function of Europa's orbital position around Jupiter . These stress fluctuations are thought to drive crevasse propagation on the moon, which can be modeled through LEFM. On Europa, diurnal components of the tidal potential (scaled over one orbital period of the moon, 3.55 days) are due to nonzero eccentricity and nonzero obliquity. Among other models, Jara-Orué and Vermeersen (2011) compute global deformation tensors from the application of Laplace transform-based normal mode theory to tidal responses. We use the analytical formulas of time-dependent tensor representing surface stresses provided by Jara-Orué and Vermeersen (2011) and due to nonzero eccentricity of 0.0094 (Wahr et al., 2009) and nonzero obliquity of 0.1 • (Bills, 2005). We neglect the contribution of relaxation modes, because the time scale of the simulation in the elastic regime is significantly smaller than the Maxwell time we adopt: 9,000 years. Further details of how the tidal stress 10.1029/2019GL084033 is computed and of the interior model are provided in Supporting Information S1. In terms of rheology, the elastic crust is sufficiently represented by Young's modulus and Poisson's ratio. We use values of 9.29 GPa and 0.33, respectively, appropriate for water ice and consistent with previous studies (e.g., Hoppa et al., 1999a, Rhoden et al., 2010. In fracture mechanics, three different cracking modes are possible, depending on the local mechanisms of deformation (Tada et al., 2000). In this study, we limit ourselves to Mode-I fracturing. This type of propagation assumes that crack flanks are moving apart from one another under stresses normal to the fracture plane. The current model is based on the calculation of the Mode-I stress intensity factor (also known as K I ) at the tip of an opening ice crevasse. Crack length, geometry, and external loads contribute to the stress intensity factor buildup. Eventually, extreme variations of these parameters lead the material to reach failure strength, called material toughness, that follows crack propagation. LEFM is usually adopted to calculate equilibrium length of cracks when this threshold is reached. Examples of LEFM application to glaciers on Earth are Hulbe et al. (2010) and Van der Veen (1998a, 1998b. These works use this technique to compute propagation of (surface and bottom) crevasses on terrestrial ice sheets, when these are loaded with longitudinal stress, overburden ice pressure, and hydrostatic pressure of water-filled crevasses.
On Europa, stress intensity factors of water-free crevasses consist of two main components: the first caused by tensile stress due to tides and the second by overburden ice pressure. If we restrict the simulation to cycloids at the surface, we can neglect the contribution of overburden pressure. Therefore, the stress intensity factor at the tip of a crevasse can be written as (Tada et al., 2000) K where is the tensile stress at the tip of the crack and l the crack length. F is an analytical function which is related to a centered crack test specimen and based on interpolation of stress intensity factor curves derived in laboratories (e.g., Tada et al., 2000; the formula is provided in Supporting Information S1). The tensile stress at the tip together with the length of the already existing crack determines whether K I reaches the material toughness of ice leading to crack propagation. If the toughness is hit and the crevasse grows in time, the increment in length l influences the stress intensity factor at the new tip. If the background stress is time dependent, at the new tip, the stress state will also be different from the previous. In the current model, we fix material toughness at 10 KPa m 1/2 , appropriate for water ice (Rist et al., 2002;Van der Veen, 1998a), and we calculate K I for advancing tips along existing cycloids.
Following Larour et al. (2004b), we adopt the displacement formulation of K I factor to compute jump-arrest propagation rates. Beside stress intensity factors, LEFM provides tools that calculate opening widths of cracks. For example, in the centered crack test specimen under stress control, the opening width at the crack center can be written as (Tada et al., 2000) = where E is the Young's modulus of the material (here, ice). V(l) is a geometrical function similar to F in equation (1) (formulas are included in Supporting Information S1).
We then derive equation (2) in time, in order to relate crevasse propagation rate to strain rate and width derivative ∕ t, which represents the opening rate. In symbols which can be rearranged in the following form: where l∕ t is the crevasse propagation rate and . the strain rate. Equation (4) implies that the propagation rate at the crack tip is a function of stress, strain rate, crack length, and opening rate. The knowledge of these four factors at the tip of the crevasse allows to estimate the crack propagation rate. This type of fracture growth is referred as stick-slip propagation with jump-arrest increment of the crack (Larour et al., 2004b). This fracturation process is characterized by extremely rapid propagation events that follows phases 10.1029/2019GL084033 of inactivity due to the conditions that are unfavorable to allow further growing. This approach was able to estimate both the nearly instantaneous cracking events that characterize propagation of rifts on Ronne Ice Shelf in Antarctica and the apparent slower growth rate that is actually measured by remote sensing (Larour et al., 2004a(Larour et al., , 2004b. There are different ways to measure the opening of a crack on ice. For example, Larour et al. (2004aLarour et al. ( , 2004b measure opening widths and opening rates of Antarctic rifts with InSAR data. Since observational constraints on surface deformation of Europa are lacking, we approximate opening width of cycloids with the displacement at the center of the crack derived from LEFM theory, that is, equation (2). Furthermore, we approximate the opening rate as follows: where . is the local strain rate. For the deformation being elastic, strain and stress solved on a plane are proportionally related by the compliance matrix, function of the Young's modulus and Poisson's ratio (Timoshenko & Goodier, 1970): where the coordinates x and y refer to zonal and meridional directions while E and are Young's modulus and Poisson's ratio. Since the elastic compliance matrix is not time dependent, stress rate and strain rate are related by the same proportion of stress and strain. Thus, time derivation of the stress tensor of Jara-Orué and Vermeersen (2011) is sufficient to compute strain rate for any location at any orbital position of Europa. However, equations (2) and (5) are first-order approximations of crack width and opening rate. In the current model, opening rates reach a magnitude of millimeters per year. On the other hand, the opening width is calculated as few meters, which is lower than the geological measurements of Galileo (Coulter et al., 2009;Kattenhorn & Hurford, 2009).
We develop a local approach that computes a series of stress intensity factors along an existing lineament, in order to calculate the growth rate of cycloids and the dormant period between following propagation events. To achieve this, we map cycloidal arcs and cusps at their exact location in the QGIS Geographic Information System software as a series of nodes separated by equidistant segments within each arc. The nodes that discretize the cycloids are obtained by using the equidistant cylindrical sphere projection of the U.S. Geological Survey global mosaic of Europa obtained from Voyager and Galileo data (Courtesy of Astrogeology Science Center, U.S.G.S, 2002). After we obtain latitude and longitude of each node through the shaping process, orientation and length of each segment are calculated with standard spherical geometry formulas (Wertz, 2001). We then compute deformation for each node of the mapped features by rotating stress tensors from the global geographic coordinate system (Jara-Orué & Vermeersen, 2011) to local referential components, using Mohr's circle (Timoshenko & Goodier, 1970): where n and t are tensile and tangential stress whereas x , y , and xy are zonal and meridional stress components, obtained from Jara-Orué and Vermeersen (2011). The angle represents the orientation of each segment, which can be computed by spherical geometry routines from the coordinates of the nodes (Wertz, 2001). Figure 1a shows one of the observed equatorial cycloids where the rotation of the stress tensor from geographical to local direction of a single node is sketched in the inset. In the figure, the rotation angle characterizes the orientation of different segments; thus, the same angle is used in equations (7) and (8).
By fixing direction of propagation and initial coordinates of the first node, the model computes stress intensity factors of advancing nodes along observed crevasses. Stick-slip propagation rates with jump-arrest increment of the crack are computed if the ice toughness is reached at the advancing node. If this is the case, the crack grows in length along the prescribed segment. Figure 1b shows a sketch of the model setup. The model calculates intensity factor, opening width, and propagation rate at N-1 advancing nodes, which are related to equidistant segments. While the background stress is constant along a single segment, different segments (at different locations) are subjected to diverse stress state. In order to determine whether the segments adjacent to the tip are contributing to the length of the open crack, our model calculates whether these are forced with tensile stress. If this is the case, we assume that the segments behind the advancing node do influence the crevasse length. This process conserves the variation of stress at different locations of Europa and the different orientation of the segments. Thus, the model repeats the analysis until the entire shape of the observed cycloid is reproduced or one of the node is not experiencing conditions that are allowing further propagation.

Results and Discussion
The cycloids observed on Europa represent an ideal case for the application of our fracture model to brittle ice. The extreme cold temperature at the surface and the consistency of the cycloid orientations with tensile stress justify our assumptions of elasticity and opening fractures. The new model is capable to reproduce the prescribed shape of four cycloids observed on the surface of the moon. We model two features at low latitudes and Delphi and Cilicia flexi at the south pole. We name the equatorial cycloids as EQ1 and EQ2. EQ2, Delphi and Cilicia are also reproduced in Rhoden et al. (2010). We discretize the four cycloids with 10,000 nodes and equidistant segments on the order of meters. Table 1 summarizes the results of our LEFM computations in terms of maximum propagation properties of the four features. As already mentioned, our model reproduces a single cycloid by bursts of propagation events, happening at each node. In the table, V indicates the maximum propagation rate at which the cycloid evolves, according to the model. This speed represents the maximum rate of a single propagation event. ∕ t is the maximum opening rate and the maximum crack width. The table includes the time required to the complete evolution of the observed cycloid and the time of inactivity. Finally, the growth rate is calculated by dividing the total length of the cycloid by the time required to complete the feature.  The growth rate is orders of magnitude smaller than the propagation rate that characterizes the bursts of activity. Propagation rates reach values of hundred of meters per seconds, while the apparent growth rate is of tens of meters per day. This is due to the fact that certain nodes require time for the stress field to rotate in order to reach the material toughness and to allow propagation. For example, Figure 2 shows a detail of the model acting on the cusp. Extreme variations in the orientation of the segments across the cusp affect the length of the open crack and the orientation of the normal stress. This means that at the cusp, the crack remains dormant. The propagation is initiated again if the stress intensity factor reaches the material toughness. Provided the observed paths, we are able to reproduce the entire geometry of EQ2, Delphi and Cilicia. The propagation of EQ1 arrests at Node 4759, and the stress field prescribed by the current tidal setup does not permit further propagation. Our model predicts that the cycloids in the southern hemisphere have a similar behavior.
Here the bursts of activity reach very high propagation rates with relative slow growth rate. On the other hand, EQ2 grows faster but with slower propagation events. This discrepancy is mainly due to the difference in the strain rate at the two locations on the surface. This is strictly prescribed by the stress forcing adopted in the model.
As already mentioned, the same behavior that our model predicts for cycloids on Europa is observed on rifts in Antarctica. Here cracks grow at an apparent rate that is slower than episodic propagation events (e.g., Banwell et al., 2017. On Earth, rifts can remain dormant for years and generally grow at rates on the order of meters per day Walker et al., 2013;Walker & Gardner, 2019). Observation from remote sensing suggests that rifts on terrestrial ice shelves experience variations in the propagation rate that often have a strong seasonal dependance Walker et al., 2015). Furthermore, several works propose that the mixture of iceberg fragments, sea ice, and snow (often called as ice melange) that fill the crack can have a substantial influence in the propagation of rifts Khazendar et al., 2009;Larour et al., 2014). The filling of crevasses on Europa should also play a substantial role in the determination of how the crack grow. Differently to what happens on Earth, Europa lacks of key observation of crack behavior and physical characterization of rifts and their filling.
However, future explorations missions will observe the moon in order to address prominent questions about the nature of the surface and interior and the prospect for material exchange between the two. For example, National Aeronautics and Space Administration's Europa Clipper (Phillips & Pappalardo, 2014) is designed to perform about 45 flybys over the nominal mission lifetime of 3.5 years. This means roughly one flyby every month, with only few opportunities for overlapping observations. The first science campaign of Europa Clipper will cover the anti-Jovian hemisphere with clusters of data overlaps during two subphases of about 6 months each (Lam et al., 2018). Therefore, overlapping observations of the same areas should be less than 6 months apart. We estimate that crevasses evolve on an average time frame of around hundreds of days. This period is roughly similar to the total observation period of a single hemisphere of the moon, designed to last less than 1 year. Average opening width would be of a few meters (value in Table 1). Therefore, we conclude that Europa Clipper should be able to detect partial development of crevasses provided that the imaging resolution of the spacecraft is of meters per pixel.
Differently to previous studies on cycloids, we model dormant periods between growth events. As consequence, our model estimates that cycloids on Europa grow at rates that are 1 to 2 orders of magnitude lower than past studies (e.g., Hurford et al., 2007, Rhoden et al., 2010. By introducing periods of inactivity, the model is capable to calculate the propagation rate at the exact location. This is a substantial difference with past models, which require longitudinal shift in order to fit some cycloids, including Delphi and Cilicia (Hurford et al., 2007;Rhoden et al., 2010). This shift was explained by secular rotation of the crust such as nonsynchronous rotation. Our model suggests that dormant period between growth events could improve fits of lineaments at their exact locations on the moon. Finally, it could also alleviate the necessity of secular motion of the crust to explain the presence of cycloids.

Conclusions
Our model extends LEFM techniques to estimate horizontal propagation scenarios by using the observed shape of cycloids on Europa. We find that crevasses evolve in a relatively short amount of time by a series of nearly instantaneous fracturing events (order of hundreds of meters per second), similarly to what is observed on rifts in Antarctica. We introduce the simulation of periods of inactivity of these cracks. This capability was not considered by previous fractures models which calculate much faster growth rates, as consequence. Beside improving fits of existing cycloids, the introduction of this critical feature could have important implication in our understanding of the rotation state of Europa. We argue that it might also lead to explain the presence of cycloids even without considering secular motion of the crust.
Furthermore, we suggest that the Europa Clipper spacecraft should be able to detect propagation events provided that the flyby repetition time over the same area is of few months apart. Our findings estimate that cracks should evolve in hundreds of days which is a roughly similar timescale compared to the coverage of the same area by the most recent trajectory baseline of the spacecraft (less than 1 year). In addition, we estimate that opening width of crevasses measures a few meters. Hence, we suggest that an imaging resolution of meters per pixel would be required to capture the scale of these events.