Frictional faulting theory is a set of equations that relate the three principal compressive stresses and that when satisfied, indicate if and what type of faults can form (Zoback, 2010). The equations are derived from the Coulomb failure envelope and Mohr circle (Anderson, 1905; Jaeger et al., 2007) and are therefore also dependent on the friction coefficient and pore pressure. The theory has been graphically illustrated in Zoback (2010) as a polygon with edges defined by the equations where points along the edges of the polygon define stress combinations associated with normal, strike-slip, and reverse or thrust faulting (Figure 2A). Reverse and thrust faulting are represented along the upper edge, which (neglecting pore pressure) is defined as S H S v = ( μ 2 + 1 + μ ) 2 , where μ is the coefficient of internal friction, 0.6 for most rocks (Byerlee, 1968). This equation simplifies to S H = 3.1 S v , when μ is 0.6 and thus plots as a horizontal line in Zoback’s polygon when S_v is known and constant. This equation can also be interpreted to state that whenever S_H is 3.1 times greater than S_v, the Mohr circle with diameter defined by these endpoints is tangential to a failure envelope defined by a μ of 0.6. Because thrust faults are observed on Mercury, (S_h, S_H) must plot as a point on the line defined by the equation above. The possibilities for the location of that point are limited due to the additional observation of the flood volcanic nature of the smooth plains units.

Tensile cracks can form in compression, parallel to the greatest principal stress (Schultz, 2019). The Griffith Criterion, which defines the stresses necessary for crack formation, allows for the relation of the maximum and minimum principal stress when cracks are opening (Hoek and Martin, 2014). Hoek (1965) found that as opening mode cracks lengthen compared to wingtip cracks, a maximum to minimum principal stress ratio of 5.88 is approached, limiting the maximum principal stress (steep blue line, Figure 2). For rocks with lower friction angles, this ratio may be even less (Hoek and Martin, 2014). An equation for the relationship between S_h and S_H when tensile cracks are observed is also derived from the Kirsch equations, which are most often used for predicting borehole stresses during drilling (Zoback, 2010, shallow blue line, Figure 2). Although the tensile fractures for which the Kirsch equations are applied are specifically produced during drilling, the stresses that produce the fractures directly at the borehole wall are often reflections of the far field stresses. The application of this second equation also narrows the window of possible stress relationships. This equation, S H = 3 S h − 2 P p   o r     S H = 3 S h .

When neglecting pore pressure (P_p), closely follows the frictional faulting equation (and edge in the Zoback polygon) for strike-slip faults. Note that by increasing the slope of this line to 5.88, one can show the relationship between stresses for the Griffith Criterion. The window, or range of stresses between the Kirsch equation and the Griffith Criterion overlaps the polygon. This window is visualized as the light blue zone overlapping the polygon in Figure 2. Most of this window overlaps with the strike-slip portion of the polygon indicating that tensile cracks can propagate in transpressive stress settings. However, some portion of this window intersects the upper edge of the polygon describing stresses associated with reverse and thrust faulting. This treatment of tensile crack propagation does not address stresses under which cracks can open in extensional stress domains.

It follows that in theory, the stress relationships described by the intersection of these two equations must closely reflect the relative stresses associated with global contraction and effusive volcanism: S_H >> S_h ≥ S_v. Dikes must be able to propagate in a compressive stress regime when the vertical stress is minimum but similar in magnitude to the minimum horizontal stress.

Dikes must be able to propagate from depth and produce enough magma to be a viable source for a large volcanic province. While Menand et al. (2010) were able to mathematically show that their equation would scale for an Earth-scale effusive event, we aimed to test if the physically derived equation for dike propagation would hold for a well-known Earth analogue site to the Smooth Plains: the Columbia River Flood Basalt and Yakima Fold Provinces (Plescia and Golombek, 1986; Watters, 1992; Crane, 2020).

The smooth plains of Mercury, and Borealis Planitia, the northern smooth plains, in particular, share important similarities with the Columbia River Basalt and Yakima Fold Provinces. The smooth plains are volcanically emplaced units with mafic, presumed basaltic, composition covering a laterally extensive region (5.59 × 10⁶ km²) in ∼0–2 km of volcanic material (Head et al., 2011; Denevi et al., 2013; Ostrach et al., 2015). These units were emplaced through multiple episodes of effusive volcanism over ∼100 Ma (Ostrach et al., 2015), a short span compared to the long-lived processes of impact cratering and thrust faulting. Northern plains deposits cover 6–7% of Mercury’s surface (Head et al., 2011) and cover the northern hemisphere in what have been inferred to be low magnesium basalts (Stockstill-Cahill et al., 2012; Weider et al., 2012; Weider et al., 2015). Thrust fault-related landforms propagate through the plain units representing post- or syn-depositional tectonics (Figure 3A). Fault propagation depths are estimated to be ∼2.5 km in Suisei Planitia, a terrain within the northern plains (Crane, 2020); however, other fault models that do not consider the effects of folding derive much greater depths (>15 km and >30 km, Peterson et al., 2020 and Egea-González et al., 2012, respectively). Comparing the emplacement time estimates for the smooth plains with the strain rate estimates for various periods in Mercury’s geologic history, we infer that emplacement occurred during the most intense periods of global contraction—the Calorian and Tolstojan periods (before 3.25 Ga, Spudis and Guest (1988)). After these periods, stratigraphic relationships with craters show that large scale volcanism ceased, and tectonism slowed substantially (Crane and Klimczak, 2017).

The Columbia River Flood Basalt Province is an extensive flood basalt province often tied to crustal extension and melting due to the mantle plume also associated with the Snake River Plain Hot Spot Track (Camp et al., 2013). However, the northern section of the basalt province is deformed through thrusting and folding of its sub-basalt, basalt, and intra-volcanic strata (Reidel et al., 2013a). Multiple effusive events covered ∼210,000 km² of Washington, Oregon, and Idaho in relatively thick (∼4 km) volcanics of largely basaltic composition (Reidel et al., 2013b). The basalts were emplaced over a geologically short period of time in multiple events between 16.8 and 5 Ma, with dikes exposed in multiple locations (Figure 2C; Camp et al., 2017). Most importantly for this study, the basalts were deforming contemporaneously with deposition in the Yakima Fold Province, a structural region of southern Washington contained within the boundaries of the Columbia River Flood Basalt Province (Figure 3B; Reidel, 1984; Kelsey et al., 2017). Faults that cut the basalts here extend to 4 km depth (Casale and Pratt, 2015; Crane and Klimczak, 2019b). In this region, uplifted, folded basalts show thinning of deposited units toward topographically elevated fold hinges (Figure 2B) and thickening in synclines and structural lows, indicating syn-volcanic thrusting and deformation (personal observations during fieldwork). Estimations of uplift from aged terraces and structural models indicate that uplift and volcanism were largely coeval, and that the intensity of both processes decreased with time (Bender et al., 2016; Staisch et al., 2018).

The Columbia River Flood Basalt Province is also a superb analogue because of the abundance of mapped dikes associated with basalt emplacement. Three main dike swarms have been described: the oldest, Steens Mountain dikes which sourced the southernmost flows, the younger, Monument dikes to the northwest of Steens Mountain, and the youngest and largest Chief Joseph dikes to the north (Barry et al., 2013). The Chief Joseph dike swarm covers northeastern Oregon and eastern Washington and was responsible for multiple basalt flows that cover the Yakima Fold Province and are deformed by its thrusts (Morriss et al., 2020). The depths to sources for the dikes are largely constrained to mid- to lower crustal sources. Davenport et al. (2017) observed a high velocity region at ∼20 km depth below a small section of the Chief Joseph dike swarm in their geophysical study of the Western Idaho Shear Zone that could be inferred to be a source. Morriss et al. (2020) estimate depths between 15 and 30 km, but state that these depths are simply the half and full thickness of the crust in the study region. Other geophysical studies constrain sources to similar depths (Hales et al., 2005; Wolff et al., 2008).

A new digitized dataset compiled from the fieldnotes of Columbia River Basalt researcher Dr. W. H. Taubeneck contains the location and description of many identified dikes. Combined with other vectorized Columbia River Basalt dikes from other dike swarms, Morriss et al. (2020) have produced a digital dataset that we used to test the physical dike propagation model of Menand et al. (2010).

We utilized a Geographic Information System and ESRI’s ArcMap to perform analyses on the dike data. All dikes and dike swarms from the dataset were combined into a single shapefile in ArcMap for analysis. The Data Management Toolbox was then used to divide each dike into 500 m long separate segments. Each segment was assigned a unique identification number and row in the new shapefile. Dividing dikes, especially very long dikes (>100 km), into segments allowed us to reflect the changing geologic conditions that might lead to varying propagation distance estimates along the length of a single dike.

Parameters were then added to this shapefile. The physically modeled equation for propagation distance relies on estimates for density of the dike magma and tensile strength and density of the dikes’ host rocks. Dike composition (usually basaltic but sometimes andesitic basalt) was recorded in the Morriss et al. combined dataset, and 2.7 g/cm³, an average density for flood basalt lavas was used (Hartley and Maclennan, 2018). Accepted values for the density of basaltic and andesitic-basaltic lavas were added to the shapefile attribute table for the dike segments. The statewide, vectorized geologic maps for Washington, Oregon, and Idaho were also imported into the GIS.

The geologic maps were used to identify the rock type through which the dikes propagated. We divided the rock mass through which the dikes propagated into a basement and bedrock layered stratigraphy. Seismic and well data for the Yakima fold Province suggest bedrock-basement contacts near 8 km (Casale and Pratt, 2015; Crane and Klimczak, 2019b). Because studies estimate crustal thicknesses between 20 and 39 km for the region (Catchings and Saltus, 1994; Das and Nolet, 1998; Davenport et al., 2017), we repeated calculations for propagation depth with three model stratigraphies to reflect maximum and minimum stratigraphic arrangements. These model stratigraphies assume that over the total distance propagated, one-tenth, one quarter, and one half of the strata are bedrock and nine-tenths, three quarters, and one half are underlying basement, respectively. We only assumed relative thicknesses for each. The equation for dike propagation distance was modified to reflect these two layers by scaling the contribution of each layer to the distance, similar to the calculation of a weighted sum.

Basement composition was assumed to be gabbro west of the Hite fault and granite east of the Hite fault, the boundary between accreted ocean crust and the North American craton (Reidel et al., 2013a). Above the gabbro or granite basement, we defined the bedrock composition as that which was recorded in the geologic map. When shallow surficial deposits were recorded as the geologic unit (i.e., alluvium or till), we instead used the host rock data recorded in the Morriss et al. (2020) dataset. The values for tensile strength and density for each of these rock types were taken from engineering sources described in Supplementary Table S1 in the Supplementary Materials. When the geologic map units were generic, we used values for rocks which most closely aligned with these descriptions. For example, values for parameters associated with shale were used when the geology was defined as “offshore marine sedimentary rock”. Tensile strength and density values for the bedrock and basement were added to the shapefile attribute table in the GIS. Once these attributes had been added, the shapefile was exported as a spreadsheet.

Dike propagation distance, d, was calculated for each dike segment. In order to calculate d, we first needed to estimate the vertical stress S_v, associated with thrust fault propagation. Based on previous estimates, thrust faults were believed to propagate upward from ∼4 km or deeper (Casale and Pratt, 2015; Kelsey et al., 2017; Crane and Klimczak, 2019b). At any depth, vertical stress can be estimated using the following formula for lithostatic stress: s v = ρ g z , where ρ is bulk density of the overlying rock, g is gravity, and z is depth (Jaeger et al., 2007). Although S_v can deviate from this value and principal stresses can rotate due to local heterogeneities such as magma chambers, some previous studies have not observed these features in the seismic data for the province and others only observe them locally (Wolff et al., 2008; Crane and Klimczak, 2019b). We therefore use the equation above to estimate S_v. For the Yakima Fold Province, ρ = 2,700 kg/m³, g = 9.81 m/s², z = 4,000 m, provide an estimate of S_v of 106 MPa. Because we know from frictional faulting theory that differential stress is estimated to be 2.1 S_v, we could predict that the differential stress at depth at the time of fault propagation in the Yakima Fold Province was 222 MPa. With the calculation of differential stress, an estimate for d could be provided for each dike segment.

The d value was added to the shapefile so that dike segments could be color coded by their propagation distance. The distances and trends in distances for dikes in the Columbia River Basalt Province were compared to previously established estimates. Similarities (discussed in section 3 Results) supported the application of the equation for propagation distance to Mercury’s Smooth Plains.

Dike propagation distance was calculated for Mercury’s smooth plains using the equation for d. Gravity on Mercury was taken to be 3.71 m/s². Tensile strength for the presumed basaltic surface composition was taken to be 14.9 MPa (Van Noort et al., 2017), although we acknowledge that depending on the degree of fracturing due to impact cratering, this value could be as low as 2 MPa (Schultz, 1993). The range in tensile strengths contribute to upper and lower bounds in our calculation of dike propagation distance.

The vertical and horizontal stresses at depth for Mercury were derived from our knowledge of frictional faulting theory and stresses associated with the topographic load of smooth plains emplacement. Topographic load was considered as this factor would have increased the vertical stresses for thrust faults to overcome and has been shown to affect the magnitude and orientation of principal stresses (Hooper et al., 2011; Roman and Jaupart, 2014). Vertical, horizontal, and shear load stresses associated with a uniform smooth plains deposit 50 km-wide and 4 km deep were calculated (Figure 4; Davis and Selvadurai, 1996; Dahm, 2000). The vertical stresses associated with this load were superimposed on the vertical lithostatic stress, and together, these stresses determined the magnitude of the maximum horizontal stress needed to overcome the yield strength of the rock mass. Some horizontal stress was also induced through loading, but the remainder was accommodated through tectonic stresses. Although we discuss global contraction as the cause for this stress, many other processes could have also produced horizontal stresses. The source of this stress is less important than its magnitude, for which a lower bound can be determined as 3.1 times the vertical stress.

The vertical stress at the depth of thrust fault initiation (2.5 km) due to a rock mass with a bulk basalt density of 2,850 kg/m³ is ∼64 MPa, and horizontal stresses associated with loading and tectonic activity were estimated to be ∼39 and ∼160 MPa, respectively. Differential stress was therefore calculated to be ∼135 MPa. Magma density was assumed to be 2,800 kg/m³ based on estimates for Earth analogue lavas and lunar basalt melts (Delano, 1990; Greeley et al., 1998). This higher density is more appropriate for Mercury, as it is believed that the surface composition of Mercury may be closer to a komatiite composition (Head et al., 2011). Applying the given parameters, upper and lower bounds for propagation distance were calculated.

Model results approximate current estimates for dike propagation distances in the Columbia River Basalts. Mid-crustal estimates from Davenport et al. (2017) and Morriss et al. (2020) range from 15 to ∼20 km while lar2gest d values in this study are ∼33 km, with most d values between 10 and 25 km. The largest values are associated with dikes that propagated through gabbro basement and meta-volcanic bedrock. Most large values (20–25 km) are generally instead associated with granite basement rock.

Trends in d values show a shallowing to the north (Figure 5). Propagation distances between 20 and 30 km were more commonly observed in the Steens and Monument dike swarms while values in the Chief Joseph dike swarm were closer to 7–10 km. Outside of the extent of the Columbia River Basalt Province, dike propagation depth variability increases with few trends in dike depth associated with direction.

Dike propagation distances for Mercury’s smooth plains were variable depending on the tensile strength for basalt used in the calculation. For intact basalt with a tensile strength of 14.9 MPa, d values ranged from 81.3 to 88.7 km. For heavily impacted basalt with a reduced tensile strength of 2 MPa, d values ranged from 10.8 to 10.9 km.

The application of an Earth analogue to this research confirmed that theoretical and physical models of dike propagation could be applied to Mercury. Dike propagation distance estimates for the Columbia River Basalt Province suggest a lower to mid-crustal source for the basalts consistent with geophysical evidence (Hales et al., 2005; Wolff et al., 2008; Davenport et al., 2017). We also observed a trend in source depths shallowing to the northern region of the province. The northward progression of decreasing source depth reflects the observed age progression of dikes younging to the north (Barry et al., 2013). These observations are consistent with the prevailing hypotheses for the origin of the basalts.

Evidence that such a large volume of basalt erupted in a short period of time and in a setting with low rates of back arc extension supports a plume origin for the basalts. A plume hypothesis holds that as the North American craton moved southwest over a mantle plume, the head and tail of the plume separated. While the tail of the plume moved along the Snake River Plain hot spot track, the head of the plume promoted melting and formation of the flood basalts (Camp and Ross, 2004). Plume models can include a delamination of the upper mantle and lower crust. The convecting plume head migrated northward and surfaceward, producing shallower crustal melt sources for the basalts over time. Segments of delaminated lower crust could have been entrained and melted within the plume head consistent with compositional data from the basalts (Camp and Hanan, 2008). Alternatively, a plume fed magma chamber with crustal contamination of craton material could produce similar compositions (Wolff et al., 2008). In either case, the high magnesium, high silica compositions of all but one flow unit imply a melted, mafic crustal source. Advancing plume heads in these models have also been connected to the development of shortening structures in front of the head (Camp and Hanan, 2008).

Other hypotheses for the flood basalt origins also align with the younging and shallowing northward pattern in the dikes. Liu and Stegman (2012) propose a model with an opening slap tear (tearing to the north) contributing to the production and northward migration of basalt. Others contend that the passive, northward propagating back-arc extension associated with the subduction of the Farallon slab would have promoted the volcanism (Carlson and Hart, 1987).

Despite an overall northward propagation pattern, the dikes indicative of the deepest source are not located near the Steens flows. Instead, they are located farther north, in the southern region of the Chief Joseph dike swarm (red points, Northeastern Oregon, Figure 5). The high d values calculated for this region are located within the Wallowa Mountains. While some studies conclude that this particular area represents the location of a plume fed and crustally contaminated magma chamber (Wolff et al., 2008), others contend that this area, characterized by at least 2 km of uplift, represents the relative zone of interaction between a plume and the granitic root of the Wallowa Mountains (Hales et al., 2005). These authors hypothesize that a plume triggered a gravitational instability, encouraging the delamination of the plutonic mountains and their subsequent uplift. In either case, these studies imply that this region may be a center of dike propagation, and so it follows that it should be characterized by the deepest source.

Our results likely reflect these important spatial patterns due to the compositions of the basement rocks and the surface rocks which formed as a consequence of these processes. If the depths had been strictly divided based on basement rock, then we would expect to only see two general dike propagation distances. Instead, the variability to the north outside of the extent of the basalt province and the deep propagation depths near the Wallowa Mountains indicate that bedrock composition also plays a role in determining the dike patterns. Smaller d values are associated with sedimentary bedrock or andesite volcanic bedrock over gabbro basement rock, due to the smaller tensile strength of sedimentary rocks and large density difference between andesite and gabbro.

Dike propagation depths for Mercury have consequences for understanding potential sources for the effusive deposits and the geologic history of the planet as a whole. When tied to compositional data, they allow us to speculate on the sources and causes of melting. Our upper bound estimates of ∼81–89 km likely correspond to some depth within Mercury’s upper mantle early in the planet’s history (Figure 6). Mantle melts would be expected to be higher in magnesium than is observed, at least in the northern plains deposits (Weider et al., 2012; Weider et al., 2015) where lower magnesium and an inferred majority plagioclase composition are observed (Namur and Charlier, 2017). This implies that for the northern plains either a thicker crust had already developed, or by analogy to the Columbia River Basalt generation, some kind of delamination had occurred to melt crustal material at mantle depths.

These depths are also consistent with depths estimated by Wilson and Head (2008). These authors considered dike propagation as function of density, crustal thickness, and the distance a dike would have to travel through the mantle before reaching the crust. They consider the stress the dike material must exert in order to keep the crack open. However, Wilson and Head present a range of solutions because this stress was unconstrained. They calculate possible dike propagation paths over a reasonable parameter space.

Our work constrains this space by calculating a narrower range of deviatoric stresses through which the dikes propagated. Our source of variability arises largely from the host rock tensile strengths. We can however learn something about mechanical stratigraphy of the lithosphere as a consequence. The lower bounds for d values (∼10 km) are not consistent with a reasonable range proposed by Wilson and Head (2008). It seems unlikely that in Mercury’s thick shell enough melt could be produced at 10 km depth to cover 40% of the planet or the vast expanse of Borealis Planitia. We take this to mean that the tensile strength of impacted basalt (2 MPa) is not indicative of a significant portion of the lithosphere. This supports the conclusions of Byrne et al. (2016) who observe nonporous (stronger) target rock and porous (weaker) target rock parameters to produce better fitting crater chronologies for larger and smaller craters, respectively. This implies that larger craters deform deeper (but still relatively shallow), stronger rock.

Global contraction is often cited as the reason that volcanism slowed (e.g., Solomon, 1978); however, finding that both processes can operate together opens up new possibilities for causal mechanisms. It is important to notice that the most intense periods of global contraction and volcanism overlap and that both processes slow contemporaneously. This implies that the shrinking lithosphere may not have been the cause for the waning volcanism. Instead, global cooling may have decreased the ability for magmagenesis and melting while simultaneously prompting contraction. The slowed rate of impact cratering would have also brought less heat to the planet, decreasing the likelihood of melting.

Viewing thrust faulting and volcanism as connected processes from a compressional stress perspective paves the way for future work studying the cessation of global contraction. Because global contraction would have peaked during the earlier periods- the time coincident with volcanism-we can assume that horizontal compressive stresses only decreased with time. Thus, by studying faults that represent a spectrum of ages postdating the smooth plains emplacement and estimating fault propagation depth, one could track the relative relaxation of stresses by calculating the deviatoric stresses associated with each fault population. It may also be possible to study pre-plains fault populations and calculate minimum deviatoric and horizontal compressive stresses for older thrusts and early periods of global contraction on Mercury.

The depths of older faults outside of the plains, but that propagated around the same time, also provide some insight into the character of diking at depth. Geometric models of faults ∼4.0 Ga place their lower tips near 30–40 km, with the elastic lithosphere extending to the same depths (Watters et al., 2002; Egea-González et al., 2012). More recent models show the present-day brittle-ductile transition (BDT) between ∼75 and 125 km (Klimczak et al., 2019). Thus, during the diking associated with smooth plains emplacement, the BDT must have ranged between these values, but closer to 40 km depth. If the maximum estimate for dike propagation distance (81–89 km) is accurate, then much of dike propagation was through the ductile regions of the lithosphere. Dikes can still propagate upward through this region, as observed in the Scandinavian Caledonides (Kjøll et al., 2019). Dikes propagating through ductile regimes have undulating contacts with host rocks, rounded shapes, and pinch and swell structures due to viscoelastic fingering and fracturing (Bertelsen et al., 2018; Kjøll et al., 2019), but importantly, dikes still migrate upwards as long as tectonic strain rates are much slower than dike propagation rates (Rivalta et al., 2015).

Because thrust fault-related landforms- and thrust faults by association-have been shown to have preferred orientations within Borealis Planitia (Crane and Klimczak, 2019a), it may be possible to map stress orientations and magnitudes through time within these plains. These stresses and their orientations and magnitudes could be linked to causes other than global contraction to explore tidal despinning and changes in orbital parameters with absolute stress values.

The methodology described here could furthermore be used to differentiate magma sources for different plains units. The differential stress value used in the source depths calculation presented in this study was derived from the fault observations in Borealis Planitia. This unit has a distinct composition from other plains units (Weider et al., 2015). At least three compositions have been inferred from spectral data (Denevi et al., 2013) and at least one of these units could represent higher degrees of partial melting and a komatiite composition- possibly implying a deeper source (Head et al., 2011). If faults in these regions are modeled, then their derived depths could be used for differential stress and d value calculations.