Finding Order in Chaos: Quantitative Predictors of Chaos Terrain Morphology on Europa.

The mechanisms for chaos terrain formation on Europa have long been a source of debate in the scientific community. There exist numerous theoretical and numerical models for chaos formation, but to date there has been a lack of quantifiable observations that can be used to constrain models and permit comparison to the outputs of these chaos models. Here, we use mapping and statistical analysis to develop a quantitative description of chaos terrain and their observed morphologies. For nine chaos features, we map every block, or region of pre-existing terrain within disrupted matrix. We demonstrate that chaos terrains follow a continuous spectrum of morphologies between two endmembers, platy and knobby. We find that any given chaos terrain's morphology can be quantified by means of the linearized exponential slope of its cumulative block area distribution. This quantitative metric provides a new diagnostic parameter in future studies of chaos terrain formation and comparison.

Introduction

Anomalous among planetary bodies in our Solar System besides Earth, Europa's surface is largely devoid of sizable (>10 km diameter) primary impact craters, with only 47 recognized in the recently released global geologic map (Leonard et al., 2019). Instead, the surface comprises numerous enigmatic features, ranging from pervasive ridges and troughs (Prockter & Patterson, 2009), to chaotic terrain resembling broken plates of preexisting terrain in a hummocky matrix (e.g., Greeley et al., 2000). The paucity of impact craters implies that Europa's average surface age is ∼60 Ma (e.g., Bierhaus et al., 2009; Schenk et al., 2003; Zahnle et al., 2003), but the exact processes and mechanisms responsible for resurfacing (erasing craters) remain poorly understood (e.g., Howell & Pappalardo, 2019). Current hypotheses addressing resurfacing on Europa range from cryovolcanism to tectonic deformation (e.g., Crawford & Stevenson, 1988; Doggett et al., 2009; Figueredo & Greeley, 2004; Greeley et al., 2000; Greenberg et al., 1999; Greenberg & Geissler, 2002; Howell & Pappalardo, 2018; Johnson et al., 2017; Kattenhorn & Prockter, 2014; Pappalardo et al., 1999), and no single proposed resurfacing model reproduces all of the current surface observations (Leonard et al., 2018).

Chaos terrain is defined as regions consisting of blocks of preexisting terrain within a hummocky matrix material. The formation of Europa's chaos has been intensely debated, with hypotheses including melt‐through, sill formation and collapse, and diapirism induced by thermal and/or chemical buoyancy (Carr et al., 1998; Collins & Nimmo, 2009; Greenberg et al., 1999; Howell et al., 2021; Michaut & Manga, 2014; O’Brien et al., 2002; Pappalardo et al., 1998; Schmidt et al., 2011). Collins and Nimmo (2009) propose a set of both hard and soft constraints on chaos formation hypotheses that are informed by observations of chaotic terrains. They note that any chaos formation mechanism must appropriately describe the matrix and block characteristics, both high and low topographies, exposure of hydrated salts and/or acids, large range in planform size (1–1,000 km in diameter), and broad global distribution. Similarly, models should consider the relationship between chaos and Europa's pits and domes (or lenticulae), preservation and destruction of preexisting terrains, and chaos growth through coalescence. Leonard et al. (2018) find that all current proposed mechanisms for chaos formation can account for only a portion of these constraints.

The ongoing uncertainty surrounding Europa's chaos terrain formation and relatively young surface age make these regions of particular interest for studies of Europa's habitability. Previous studies hypothesize that chaos terrain is an expression of the subsurface ocean interacting with the surface, whether exchange is occurring through the solid, frozen transport of icy materials (e.g., Howell & Pappalardo, 2018; Johnson et al., 2017; Kattenhorn & Prockter, 2014), or through the formation, migration, and eruption of liquid water and brines (e.g., Crawford & Stevenson, 1988; Greenberg et al., 1999; Schmidt et al., 2011). The exchange of material between the surface and subsurface of Europa may be critical to sustaining the chemical disequilibria that could allow the emergence or persistence of habitable environments and potentially life (Hand et al., 2007, 2009; Vance et al., 2016). Thus, determining how chaos terrain forms is an important key to assessing the habitability potential of Europa's icy shell and subsurface ocean. In order to constrain formation mechanisms, existing studies note observations of chaos morphology (e.g., Collins & Nimmo, 2009; Leonard et al., 2018; Skjetne et al., 2021), the geographical distribution/trends (Noviello et al., 2019), or focus on one chaos terrain (e.g., Figueredo et al., 2002; Spaun et al., 1998), but quantitative morphological observations that can be used to directly compare individual terrains remain elusive.

Existing chaos models tend to cite a few prominent examples of chaos terrain that are consistent with a modeled formation mechanism without considering the full array of chaos morphologies. The two prime chaos terrains used as type examples are Conamara Chaos (TH_cc in Figure 1) and Murius Chaos (LH_E in Figure 1). However, we hypothesize that chaos on Europa occupies a morphological spectrum, likely resulting from variations in the same formation mechanism, and that the previously mentioned oft‐used examples of chaos are more akin to endmembers on this spectrum. Thus, we would expect quantitative metrics of terrain morphology to transition smooth across a range of values as morphologies transition from platy to knobby. This hypothesis would be challenged if the result was that these quantitative metrics are uncorrelated to the observed morphologies, or by distinct grouping of end‐member terrain types without intermediate values that span them. In this work, we provide a means to quantify chaos morphologies through geologic mapping and investigation of the size‐frequency distribution of blocks within chaos terrains, facilitating comparisons. The result is a straightforward, quantitative metric permitting chaos morphologies to serve as a diagnostic tool for modeling efforts.

Figure 1

(a) Europa global mosaic with chaos terrains (black outlines) used in this study. The TH and LH labels denote the chaos terrains that are in the Trailing and Leading Hemispheres, respectively. (b) TH_cc (Conamara Chaos) with outline of chaos terrain borders mapped (left) and chaos blocks mapped (right). See Figure S1 in Supporting information S1 for images of every other mapped chaos terrain and Table S1 in Supporting Information S1 for information on all images. (c) Spectrum of representative examples of each chaos terrain mapped in this study, ordered from most platy (left, 1) to most knobby (right, 9).

Geologic Setting and Approach

Europa's chaos terrain generally consists of high relative brightness blocks of crustal material, ranging from 100 – 101 km wide that reside within a lower relative brightness hummocky matrix (e.g., Doggett et al., 2009; Greeley et al., 2000; Kattenhorn & Hurford, 2009; Leonard et al., 2020). Applying the terms defined by Greeley et al. (2000), chaos morphology is characterized as “platy” or “knobby” (Figure 1). Platy chaos is defined as a unit that consists of clearly defined “plates” of pre‐existing terrain, where the texture and morphology of the pre‐existing terrain is clearly visible due to scale or lack of disruption. Knobby chaos is defined as a unit that consists of smaller blocks, or “knobs” of pre‐existing terrain where the pre‐existing morphologies are not apparent. In this study, we identify a chaos block as a piece of terrain within the broader margins of a chaos unit that has a different texture and topographical expression (standing higher or lower based on shading and shadows) than the surrounding hummocky matrix (Figure 1b). Both platy and knobby chaos can have sharp boundaries with the surrounding terrain.

We determine the effective diameter of blocks by employing size‐frequency analyses on the blocks mapped within individual chaos units. Generally, size‐frequency distributions (SFDs) are used to relate the magnitude of an event or object type with its frequency of occurrence. Here, we investigate the frequency of occurrence of chaos block sizes within individual chaos units. This SFD reveals clues about the formation of the terrain, and allows for quantitative comparisons across different types of chaos terrains. Size‐frequency distributions are most well known in planetary science application for investigating craters on a planetary surface for use in absolute and relative dating (e.g., Crater Analysis Techniques Working Group, 1979; Zahnle et al., 2003). More broadly, SFDs can be found in a wide range of applications including characterization of terrestrial rockfalls (Ruiz‐Carulla et al., 2018), fragmented ice (Toyota et al., 2011), and Mars landing site rocks (Golombek et al., 2003; Golombek and Rapp, 1997).

In studies of rockfalls on Earth, block size distributions are employed as an indicator of fragmentation degree. There are several parameters that can affect fragmentation in rockfalls including the presence of discontinuities, impact angle, impact energy (height), and rigidity of the ground. By measuring the size distribution and some of the variables that affect fragmentation, Ruiz‐Carulla et al. (2018) find that the rockfall block size distributions are strongly correlated to the height of the fall. Likewise, iceberg field SFDs are predicted to relate to the fragmentation process, (i.e., the fracturing process that separated the chunk of ice from the glacier) (Åström et al., 2021) and attempts have been made to apply this concept to Europan chaos (Walker & Schmidt, 2014). However, by studying the block size distribution of different chaos terrains on Europa, we may be able to provide observational constraints on formation mechanism and hypothesize what factors control the size distributions.

Methods

Chaos Characterization

In this study, we focus on regions where chaos terrains have previously been identified (e.g., Figueredo & Greeley, 2004), where chaos terrains are fully contained within regional resolution imaging (200–300 m/pixel), and where the total extent of the chaos terrains are at least 50 km in diameter. This size constraint is chosen to be large enough to ensure that numerous blocks will be present within the chaos terrain. Moreover, it is essential that the borders of the chaos terrain are fully visible in the regional resolution imaging, as to not potentially skew the resulting SFD; we do not use partially‐imaged chaos terrains as they may not be representative of the terrain as a whole. Additionally, we choose a regional resolution of 200–300 m/pixel in order to fully capture blocks down to 2 km in its shortest dimension. We identify chaos terrains of a similar latitude (0–30° North) in order to avoid the hypothesized potential influence of latitudinal variations in tidal deformation and heat production (e.g., Nimmo et al., 2007). Applying these constraints, we perform mapping and measurements for nine chaos terrains (Figure 1a, Table S1 in Supporting Information S1). We map the borders of each chaos terrain and each block within a chaos terrain greater than 2 km in its smallest dimension (see Figure 1b). By accomplishing this mapping using ArcGIS, we can then determine the area of each block in each mapped chaos terrain.

After mapping, we order the chaos terrains by morphology, from most‐platy to most‐knobby (Figure 1c). For each terrain, we interpret the morphology and assign a morphological rank from 1 to 9, where 1 is the most‐platy and 9 is the most‐knobby, creating a spectrum of morphologies. This ranking is based on subjective mapping observations alone and was performed by two of the authors (A. Mills and E. J. Leonard). The respective rankings were then compared and a final ranking agreed on by both authors. Importantly, when directly compared, chaos morphologies occupy a spectrum of platy‐to‐knobby behavior where the two commonly cited examples of chaos terrain, Conamara Chaos (TH_cc) and Murias Chaos (LH_E) are close to end‐members on the spectrum (ranked 1 and 9 respectively). This spectrum‐like behavior implies that creating an artificial divide into two morphological groups (e.g., Figueredo & Greeley, 2004) is likely fraught with subjectivity and should be avoided when mapping.

Cumulative Size‐Frequency Distributions

In order to objectively and quantitatively describe chaos terrains and their constituent blocks, we derive a cumulative distribution for block areas (Figure 2). We first calculate a characteristic block length scale, x, which we define as the square root of the block area. We then generate cumulative SFDs by relating the cumulative area of blocks of a given characteristic length scale to the characteristic length scale of any given block, similar to an R‐plot used in cratering statistics (Crater Analysis Techniques Working Group, 1979).

Figure 2

Cumulative block area size‐frequency distribution for 9 chaos terrains mapped in this work. The data points for each chaos terrain are connected by solid straight lines in order to highlight the trends. The dotted line for each chaos terrain is the linear fit of the form given in Equation 2.

The cumulative fractional area of blocks with a length scale larger than x is represented by:

Here, is the sum of the area of every block with a characteristic length scale greater than , is the total area of the given chaos terrain, and m is a constant. represents the total cumulative area that blocks of all sizes (as goes to zero) occupy in a given chaos terrain. We normalize the fractional area of blocks greater than any given characteristic length scale () by the total chaos terrain area (), because larger chaos units have an opportunity to contain more blocks (Figure 2). For easier comparison among distributions, Equation 1 can be linearized as: where . In Figure 2, we plot the natural log of the fractional cumulative area of blocks (left side of Equation 2) versus the characteristic length scale of the blocks () for each chaos terrain mapped, along with the associated least‐squares fit. The results are described in Section 3.

Results

The data displayed in Figure 2 are best fit by exponential functions with varying constants (and of Equation 2), with residuals varying from 0.95 to 0.99 (Table 1). The slopes of the linearized exponential best‐fit lines range from about −0.2 to −1.1 in units of km−1. Likewise, the y‐intercepts of the best fit lines range from about −1.5 to −0.2, with one outlier of ∼1. These y‐intercepts represent the natural log of the percentage of area occupied in each chaos terrain by blocks of any size (down to zero), or , by definition. For example, the TH_bb has a y‐intercept of −1.0254, meaning the percentage of chaos area occupied by blocks of all sizes (as goes to zero), is exp(−1.0254) or ∼36%.

Table 1

Summary of SFD Results

Chaos terrain (Figure 2)	Chaos area total, A_C(km2)	Block area, A_x>2(km2)	Slope, m(std. error)	y‐intercept, b(std. error)	Fractional area taken up by blocks with x > 0, A_b	R 2
LH_A	12,242	4,426	−0.2712 (0.0037)	−0.3853 (0.0157)	68%	0.96
LH_B	2,765	909	−0.5129 (0.0105)	−0.2374 (0.0317)	79%	0.96
LH_C	2,945	1,582	−0.2140 (0.0020)	−0.2015 (0.0086)	82%	0.99
LH_D	1,959	667	−1.1065 (0.0195)	1.0694 (0.0536)	291%	0.98
LH_E	7,682	865	−0.7445 (0.0192)	−0.8715 (0.0551)	42%	0.95
TH_aa	13,186	2,835	−0.4377 (0.0051)	−0.7785 (0.0173)	46%	0.97
TH_bb	22,570	4,854	−0.2995 (0.0024)	−1.0254 (0.0092)	36%	0.98
TH_cc/Con. Chaos	8,997	3,393	−0.2653 (0.0039)	−0.3997 (0.0172)	67%	0.97
TH_dd	7,979	798	−0.4366 (0.0137)	−1.5420 (0.0477)	21%	0.95

In order to investigate the relationships among measured parameters, we perform a cross‐correlation analysis between the measured and derived properties of each chaos region (Figure 3). Specifically, we look for relationships among: (a) the total area of the chaos region (, (b) the morphological rank that resulted from ordering the chaos terrains based on the subjective prevalence of visible platy and knobby characteristics (Figure 1c), (c) the cumulative area taken up by blocks of  > 2 (), (d) the fraction of the area taken up by blocks of  > 2 (, which we define as ), (e) the exponential slope of the SFD (), (f) the exponential intercept (), (g) chaos center longitude (east positive), and (h) chaos center latitude (0–30 N), as defined in the study parameters.

Figure 3

Cross‐correlation matrix with scatter plots showing the relationship among mean parameter values, for each chaos region. The Spearman's rank correlation coefficient, , and Type I error rate, , are noted beneath each correlation plot (see Section 3.2). Significant correlations () are highlighted in red. Note, the and values cited for D5 do not include LH_D, an outlier.

For the correlation analysis, we employ a Spearman's rank correlation (Hauke & Kossowski, 2011), which tests for any monotonic relationship between two variables. While the Pearson's moment product is more common in statistical analyses of planetary datasets, the Pearson's moment product tests for a linear relationship between two variables that may not be suited to the exponential, power‐law, and other complex relationships among the variables studied. The Spearman's rank correlation, however, does not require that correlated variables have linear relationships. When increasing the value of one variable of interest, this approach yields a correlation coefficient, r, of 1 when always resulting in the increase of another variable, of r = −1 when always resulting in a decrease, and of r = 0 when equally likely to result in an increase or a decrease (Hauke & Kossowski, 2011). We define a significant correlation as exhibiting ≥95% probability that the relationship between two variables is not the result of random sampling of two unrelated variables, as characterized by the p‐value (p ≤ 0.05).

Figure 3 illustrates the cross‐correlation analysis, where statistically significant correlations are indicated in red. The relationships that were found to be statistically significant are: (a) exponential slope of the SFDs () and morphological rank; (b) block area () and total chaos area (); (c) exponential intercept () and fractional block area (); (d) total chaos area () and geographic location; and (e) exponential intercept () and geographic location (Figure 3). In the next Sections (4.1–4.4) we discuss the implications of these correlations.

Discussion and Perspectives

Exponential Fit

As noted in Section 3, we find that the cumulative areas of the chaos blocks are best fit by exponential functions (Figure 2). As described by Golombek and Rapp (1997), exponential functions for describing block (or rock) abundance are better suited than a power‐law distribution, considering the behavior of these functions at small sizes (e.g., Turcotte, 1986). Likewise, exponential functions are also found to be a better fit to cumulative fracture length distributions (e.g., Cowie et al., 1993; Howell et al., 2016). The exponential fit fails at small and large block sizes (Figure 2), which is expected because of likely undercounting of small blocks and the naturally occurring small numbers of large blocks.

The y‐intercept of the fits () characterizes the behavior of the uncountable and unmappable small blocks (i.e., as goes to zero), representing the percent of cumulative area that blocks occupy for all length‐scales (). These vary between −1.54 and −0.22 (Table 1). One chaos terrain (Chaos D) exhibits a non‐physical value greater than 1. Likely, this occurs because the current functional form does not accurately describe the cumulative block area for blocks of x << 2 km (see Section 4.3 for further discussion on this point).

Slope and Morphological Rank (Correlation)

The slopes of the linearized exponential SFDs (Figure 3, A4) show statistically significant correlation to the more subjective morphological rank of the chaos terrains. The correlation is such that a chaos terrain with more negative exponential slope has a higher morphological rank. Based on how we defined the morphological rank (Figure 1c), a higher morphological rank indicates a more‐knobby morphology, whereas a lower morphological rank indicates a more platy morphology. Thus, a steeper (more negative) exponential slope would indicate that the chaos terrain has a knobby morphology where as a shallower (less negative) exponential slope would indicate that the chaos terrain has a platy morphology. This correlation allows for a quantitative description of the morphology of any chaos terrain. Therefore, it is also possible to predict the morphology of a chaos terrain given a cumulative block area distribution, potentially providing a useful test to numerical models of chaos formation.

Interestingly, the SFD slopes (, Table 1) fall across a bounded and continuous range of values with no clear changes in behavior, even as the chaos morphologies change. This suggests that each chaos terrain studied may be governed by the same principles of formation, and potentially the same formation mechanism. If this is the case, the varying SFD slopes could be indicative of the varying conditions under which the same process has acted to form chaos. By analogy to previous work on icebergs and rockfalls (e.g., Åström et al., 2021; Ruiz‐Carulla et al., 2018) the correlation between the slope of the block size distribution and the morphology may be indicative of the amount of energy partitioned into fracturing the terrain when the chaos feature formed. For example, knobbier chaos regions tend to contain a larger fractional area of smaller blocks compared to the platy chaos terrains; as the ratio of small blocks to large blocks increases, more block‐bounding fractures are required to break up the blocks, increasing the energy required to form the unit.

The correlation between SFD slope and morphology can also be used as a tool when creating geologic maps of Europa's surface. Instead of mappers creating an artificial divide on a seemingly continuous morphological spectrum, the chaos SFD could be used. For example, for chaos terrains large enough to display mappable blocks (which depends on the image resolution and scale of the map) blocks should be mapped and the slope of the SFD reported in place of subjectively assigning qualitative unit descriptions. However, at image resolutions or mapping scales that do not permit blocks within chaos regions to be mapped, a qualitative description remains appropriate.

Total Block Area and Chaos Area (Correlation)

Another significant correlation is that between the total block area (for blocks with a characteristic length >2 km, ) and total chaos area (Figure 3, B2). As a chaos terrain area gets larger, the total area of the blocks within that chaos terrain also gets larger. To first order, if all chaos terrain was formed by the same process that involved breaking blocks of pre‐existing terrain, we would expect larger chaos terrains to have more blocks—and therefore more total block area—than smaller chaos terrains. However, observations of chaos reveal that the story is more complicated than a process that simply breaks apart the surface: blocks are rotated, partially entrained in a matrix‐like material, or completely missing (assumed to be submerged in or converted to matrix material). We predict that the process(es) forming chaos act in a way that similarly breaks apart, rotates, and engulfs/degrades blocks across all chaos terrains, proportionally with chaos size. Thus, this correlation may lend credence to the hypothesis that all chaos terrains are formed by the same process acting at different scales.

Block area appears to increase approximately linearly with total chaos area (Figure 3, B2). This correlation is consistent with the hypothesis that larger chaos terrains are formed by combination or merging of many smaller chaos terrains. Analysis of future higher‐resolution imagery of Europa's surface, when smaller chaos terrains cumulative block area is able to be determined, could support this idea further if it is found that adding two SFDs for small chaos together can produce a collective distribution similar to that of a different chaos terrain of a similar large size.

In direct contrast, block area fraction (, defined as ) is not found to be significantly correlated with block area () or chaos area () (Figure 3, B3 and C3). Because block area and chaos area are correlated, the expectation is that block area fraction would also be correlated with chaos area and block area. The fact that block area fraction does not show a significant correlation with chaos area or total block area could indicate that the relationship between total chaos area and block area is complex, and therefore not indicative of larger chaos terrains being formed by merging smaller chaos terrains. Resolving this conflict is therefore an important target of future work that may especially benefit from new observations of chaos terrains across different regions and at higher resolution.

We also find that (the y‐intercept of Equation 2) is correlated with block area fraction (Figure 3, D5). As the block area fraction () increases, the exponential intercept () also increases. This correlation is interesting because the exponential intercept is related to the cumulative total block area as goes to zero (, by definition) and the block area fraction (, defined as ) is the total block area of blocks with a characteristic length >2 km. For LH_D, the chaos that has a b > 1, it could be the case that the single exponential fit break down at some point for blocks with  < 2 km, indicating that some other process is dominating smaller block sizes for this chaos unit. In the Galileo E12 high‐resolution (∼6 m/px) sequence across Conamara Chaos, it appears that chaos blocks do have some minimum typical size (i.e., there is not a continuum of block sizes down to the limit of resolution), but what this minimum size is requires further investigation and high‐resolution data across many chaos regions. In general, because block area fraction is correlated to b, this indicates that the behavior of smaller blocks is still well represented by a single exponential function.

Other Correlations

The latitude and longitude of each chaos terrain unit are also significantly correlated to the slope of the cumulative block area distribution (Figure 3, F6 and F7). It has been speculated that chaos distribution and morphology might vary with latitude and/or longitude (e.g., Doggett et al., 2009). However, because we only sample two longitude bands (see Figure 1A and Table 1), these correlations are most likely an artifact of sampling bias. Future mapping of chaos terrains spanning greater latitudinal distribution may allow the latitudinal relationship to be fully investigated, while new and more extensive imaging, expected from the Europa Clipper mission, is required to permit mapping of chaos terrains at different longitudes.

Future Work

The Europa Clipper spacecraft is expected to arrive at Europa in 2031, and will supply a wealth of new imaging data of Europa (Howell & Pappalardo, 2020). In particular, the EIS wide angle camera and narrow angle camera will build up coverage over ∼50 flybys resulting in >80% of the surface imaged at better than 100 m/pixel (Turtle et al., 2019). This data set will allow for a global assessment of chaos block size distributions, releasing existing limitations on latitude, longitude, and chaos size. Additionally, the Juno spacecraft will flyby Europa in Fall 2022 and image a part of Europa which could potentially enable new work.

In the meantime, chaos block size distributions are a quantifiable observable that can be used to test chaos formation models. Models of chaos formation could explore mechanical controls on chaos morphology, identifying the key physical parameters that control SFD slope in order to provide better insight into the formation mechanism, mechanical controls, or timescales. Understanding these aspects, and others, of chaos formation could indicate the presence of water or chemical gradients in Europa's ice shell. This work thus opens up a new quantitative tool that can help us explore the relationship between geology and habitability.

We find that the slope of a cumulative block area distribution is well‐correlated with the morphological rank of a given chaos terrain and that this chaos morphology occupies a spectrum, indicating that variation within a single formation process (e.g., energy input into fracturing the surface) is the likely driver of the spectrum of chaos morphologies. This quantification of chaos morphology now allows for deeper and more direct comparisons of chaos terrains. It is not yet clear whether chaos formation involves liquid water and/or communication with the subsurface ocean, but we can now move forward with a quantifiable benchmark to differentiate among chaos formation models.

Supporting Information S1

Click here for additional data file.