No Unique Scaling Law for Igneous Dikes

In linear elastic fracture mechanics (LEFM), veins, dikes, and sills grow in length when the stress intensity factor KI ${K}_{I}$ at the tip reaches a critical value: the host rock fracture toughness KIc ${K}_{Ic}$ . This criterion is applied broadly in LEFM models for crack growth and it is often assumed that the pressure inside the crack is uniform. When applied to intrusion length versus thickness scaling, a significant issue arises in that derived KIC=300to3000MPam ${K}_{IC}=300\,\mathrm{t}\mathrm{o}\,3000\,\mathrm{M}\mathrm{P}\mathrm{a}\,\sqrt{m}$ , which is about 100–1,000 times that of measured KIc ${K}_{Ic}$ values for rocks at upper crustal depths. The same scaling relationships applied to comparatively short mineral vein data gives KIc<10MPam ${K}_{Ic}< 10\,\mathrm{M}\mathrm{P}\mathrm{a}\,\sqrt{m}$ , approaching the expected range. Here we propose that intrusions preserve non‐equilibrated pressures as cracks controlled by kinetics, and therefore cannot be treated in continuum with fracture‐controlled constant pressure (equilibrium) structures such as veins, or many types of scaled analogue model. Early stages of dike growth (inflation) give rise to increasing length and thickness, but magma pressure gradients within intrusions may serve to drive late‐stage lengthening at the expense of maximum thickness (relaxation). For cracks in 2D, we find that intrusion scaling in non‐equilibrium growth is controlled by the magma injection rate and initial dike scaling, effective (2D) host rock modulus, magma viscosity and cooling rate, which are different for all individual intrusions and sets of intrusions. A solidified intrusion can therefore achieve its final dimensions via many routes, with relaxation acting as a potentially significant factor, hence there is no unique scaling law for dike intrusions.

GILL ET AL. 10.1029/2022JB024120 2 of 13 including material fracture toughness (Scholz, 2010); cracks would become unstable under constant stress loading, therefore implying growth under constant displacement boundary conditions (Segall, 1984). Understanding scaling relationships therefore has significant implications for the mechanics of intrusions and other opening mode fractures.
Opening displacement (thickness) versus length ( -) data for dikes (and veins, sills, etc., but here we focus on dikes) are generally interpreted using a linear elastic 2D pressurized crack model. The model assumes mechanical equilibrium, such that the stress intensity, , at the tip of the dike is equal to the mode I fracture toughness of the country rock, (the ability of a material containing a crack to resist fracture). Magma flows within a conduit down a pressure gradient, so a static (equilibrium) condition necessarily requires that the magma pressure, , is uniform within the dike, as shown in Figure 1b. In reality, the ability of a magma to flow to relieve overpressure and achieve equilibrium will be directly dependent on magma viscosity, for which there is a significant range in nature (McLeod & Tait, 1999), rising sharply toward solidus temperatures. For our purposes here, in a 2D case, we will consider the length ( ) as the maximum horizontal dimension of the dike (as opposed to the dike height; the maximum vertical dimension), and the thickness ( ) as the minimum horizontal dimension of the dike ( Figure 1b). The failure condition at constant pressure is The maximum thickness would be at the center of this 2D dike and is given by where ′ = 1− 2 is the (2D) plane strain modulus of the country rock, is the Young's modulus and is the Poisson's ratio. Combining Equations 1 and 2 gives the classic relationship between dike thickness and length (Olson, 2003) where the constant of proportionality is = √ 8 ′ . Figure 1a suggests that 10 −2 < < 1 √ encompasses the range of dikes observed in the field.
Measured thickness to length ratios are generally consistent with reasonable magma excess pressure estimates using Equation 2, in the range of 1-10 MPa (Rubin, 1995), but the large areas over which that pressure operates in a constant pressure model results in extremely large stress intensity at the tip, which then requires excessively large fracture toughness to stabilize the crack (Figure 1a). It is widely acknowledged that this model-predicted value for fracture toughness is much larger than the expected values for the host rock (Cruden et al., 2017;Rivalta et al., 2015), with model results typically in the region of = 300 − 3000 √ on Earth (Schultz et al., 2008), with estimates up to = 15, 000 √ on Mars (Rivas-Dorado et al., 2021); these predictions are compared to ∼1 √ measured in the laboratory (Atkinson, 1984). To appreciate how large these calculated values are, the fracture toughness of all classes of material are shown in Figure 1c, where ≈ 10 −5 √ for technical ceramics including glass, and ≈ 10 −4 √ for building materials such as concrete and brick, with the very highest values of ≈ 10 −3 √ for high-performance structural materials such as metal alloys. Predicted dike-model values of are 2-3 orders of magnitude above the expected and measured range for rocks, and significantly above the toughest known materials, such as maraging steel (175 √ ) and titanium alloys (up to 107 √ ). Measured for upper crustal rocks (0-5 km) ranges from about 0.5 − 3 √ (Stoeckhert et al., 2015), hence the equilibrium model of Equation 1 cannot explain what is physically observed using realistic material parameters. The high fracture toughness of metals shown in Figure 1c is due to their plasticity, and it could be argued that rock plasticity and/or increasing depth/temperature should increase (Balme et al., 2004;Heimpel & Olson, 1994;Stoeckhert et al., 2015), but this is not enough to span the expectation gap, especially at the shallow crustal emplacement depths of the dikes plotted in Figure 1a. Delaney et al. (1986) proposed that the discrepancy between measured and calculated fracture toughness values relates to the difference in size of fracture process zone between the sample-scale (small) and the field-scale (large). The process zone should scale with crack length, hence longer fractures should expend greater energy in . Graben data (Elysium data form Rivas-Dorado et al., 2021; Tharsis data from Mège et al., 2003) refers to dike dimensions, based on calculations using graben widths (Rivas-Dorado et al., 2021). Sudan (Babiker & Gudmundsson, 2004), Ethiopia (Schultz et al., 2008), and Shiprock (Delaney & Pollard, 1981) dike data and Culpeper and Florence Lake vein data (Vermilye & Scholz, 1995) were extracted from syntheses presented by Olson (2003) and Cruden et al. (2017). Deccan dike data is from Ray et al. (2007) and Karoo dike data is from (Coetzee and Kisters (2017). Sandstone (SST) dike data from Vétel and Cartwright (2010): PGIC, Panoche Giant Intrusion Complex (California, USA). (b) Pressure and dike thickness profiles for toughness-controlled (upper image) and kinetic-controlled (lower image) models, showing the mode I stress intensity at the dike tips, , in both cases. (c) Materials Selection Chart (adapted from Ashby, 2009) showing versus for a range of materials. Note the position of values predicted from equilibrium-based intrusion scaling relationship models relative to the position occupied by natural rocks. Full explanation is given in Section 2. 4 of 13 required to form a process zone ahead of a hydrous vein than for a dike, at a given length-scale, despite dikes and veins being hosted in similar rock types. The GPa-scale values calculated for are for already-long intrusions, with the stress intensity proportional to the crack length (Equation 1), hence longer intrusions should be easier to grow than short intrusions; a toughness-controlled equilibrium model effectively predicts that it is impossible to grow a short intrusion, since ≪ at shorter length scales.
Any alternative model for dike growth must therefore predict a stress intensity at the dike tip that is within a realistic fracture toughness range for rocks in the upper crust-on the order of ≈ 1 √ -which is the purpose of this paper. However, the existing toughness-controlled growth model would appear to be inappropriate. As an illustration, taking a typical dike from Figure 1 with = 6 and = 1 , the equilibrium model predicts that a host rock with ≈ 1000 √ is required to sustain this dike. If we keep the magma volume (area in the 2D case) constant, and reduce the fracture toughness to ≈ 1 √ then Equation 3 predicts that this dike would have dimensions of = 6 and = 100 at equilibrium. This shape is never likely to be achieved of course, since magma flow would cease due to solidification, and the final dike would be one that is frozen into a non-equilibrium state.
Here we revisit the assumptions of dike scaling laws, reapplying principles of kinetic (viscous) and fracture controls on 2D crack growth. In particular, we assess whether dikes and veins occupy the same − scaling continuum, and examine the fundamental controls on their growth and preservation in the rock record that are manifest in such scaling plots. A range of realistic conditions for individual dike systems are modeled and show that no two systems are likely to follow the same scaling trajectory.

Kinetic-Dominated Versus Toughness-Dominated Growth
Here we will consider the growth of cracks in 2D in terms of the evolution of length and thickness, in which the length represents the maximum horizontal dimension, and the thickness is the minimum horizontal dimension (e.g., Figure 1b). Equation 1 assumes the dike is in an equilibrium state when it solidifies, that is, the fluid (magma) has had time to redistribute itself within the fracture to remove all pressure gradients. This toughness-dominated assumption is reasonable for low viscosity fluids in small fractures, since the short distances involved mean that fluid pressure can equilibrate quickly. However, as noted above, non-equilibrated fracture geometry results in the prediction of unphysically high fracture toughness, so it is necessary to look at alternative explanations. The principal variable in linear elasticity that has influence on , and that can be changed, is the excess magma pressure distribution, ( ) , where − 2 ≤ ≤ 2 is the lateral distance from the center of the dike. The exact pressure profile in the dike is represented analytically by a series expansion (Spence & Sharp, 1985) but can be illustrated in simpler terms using the approximate solution of Spence and Turcotte (1985). They introduced a linear variation in the pressure such that ( ) = + Δ | 2 | , where is the pressure at the center of the dike and + Δ is the pressure at the tip. A value of Δ = − 2 removed the stress intensity at the tip entirely, that is, = 0 , resulting in a negative tip pressure of − 0.57 , as shown in Supporting Information S1. This is illustrated in the lower diagram in Figure 1b. For a pressure that continually decreases away from the center, which is consistent with magma flowing toward the dike extremities, the excess pressure at the tip must always be negative for = 0 (note, this is only a negative excess pressure, and once lithostatic pressure is included, the total pressure is still positive and therefore compressive). This alternate extreme is referred to as kinetic-dominated behavior, whereby a dike propagates in a non-equilibrium state determined by the rate at which magma is emplaced and redistributed within the fracture. It assumes that the fracture toughness is negligible compared to the large forces involved in dike propagation, such that it can be assumed that is effectively zero. The primary assumption behind this model is that the crack tip must remain magma-filled, whereby any (low pressure) cavity that developed would quickly be filled or closed due to the large pressure difference between the magma or host rock and the cavity (Rubin, 1995). This model was first employed by Delaney and Pollard (1981), and it is generally accepted that the precise tip conditions in this respect are not that important (Rubin, 1995); for example, the existence of tip cavities that are small compared to the dike length do not change the = 0 assumption. The kinetic-dominated versus toughness-dominated argument has been discussed in the literature for some time (see Rivalta et al., 2015). It is therefore useful to quantify the predicted − scaling response where these regimes apply. To do this we utilize the simple 2D analytical approximation of Spence and Turcotte (1985), which allows for both finite toughness and finite viscosity, to model the growth of a (2D) dike with linearly increas-5 of 13 ing volume (area) 2 = , where is the injection rate (in units 2 −1 ) as a function of time. Firstly, it is of particular interest to note that this allows us to obtain the same scaling relationship as Equation 3 (see Supporting Information S1) but with a different constant of proportionality, where ( ) is a dimensionless scaling function which is a function of the dimensionless scaling parameter is the toughness-dominated length scale, and = ( ′ ) 1 2 is the kinetic-dominated length scale. In the latter, is the constant growth rate ( 2 −1 ) and is magma viscosity ( ). The parameter is the key measure of the balance between toughness-controlled ( 1 ) and kinetic-controlled ( 1 ) growth. A plot of ( ) versus is shown as the blue line in Figure 2. It shows that = 1 where toughness dominates, as expected from Equation 3, and that 1 where kinetic effects dominate. The red line shows the predictions for the purely kinetic regime ( = 0 ) for which in contrast to the toughness-dominated prediction of Equation 3. It is clear that systems are kinetics-dominated for 0.2 and toughness-dominated for 0.4 where = 1 . There is a small transition region in between, but it appears this is not significant enough to warrant a combined model; that is, it is sufficient to use a toughness-dominated or a kinetics-dominated model. It is useful now to estimate where a particular system sits on this continuum, particularly when reflecting on the scaling relationships between laboratory models and natural fracture systems.
For dikes, if we assume = 10 −2 − 1 2 ∕ , = 10 2 − 10 8 , = 10 6 √ and ′ = 1 − 10 , then we get 10 −4 < < 10 −1 (Figure 2). Spence and Turcotte (1985) used lower viscosities and estimated ≈ 5 × 10 −3 , but in any case, it is clear from Figure 2 that dikes are very strongly dominated by kinetics. Parameter can have a wide range of values, depending on the specific conditions under which a dike was emplaced, and predicts, therefore, a wide range of observed dike aspect ratios, consistent with the wide GILL ET AL.
10.1029/2022JB024120 6 of 13 scatter in the observed data. This model suggests that rapid emplacement (large ) of a viscous magma (large ) into a compliant host (low ′ ) leads to the growth of a relatively short and thick dike (small and large ). The chosen viscosity range is high for a basaltic magma (typically taken as = 10 2 ), but in line with phenocryst-rich andesite or rhyolite magmas (Takeuchi, 2011 and references therein). Viscosity is a strong function of temperature, hence the viscosity of even basaltic magmas will approach such a high value as they approach solidus temperatures; a condition that becomes more likely toward the periphery of an intrusive system. In this model, low pressure gradients and slower plug-flows would reduce the effective channel width of the conduit, consistent with a higher viscosity.
In vein systems formed by hydrofracture (for purposes of comparison to Figure 1a, we are referring exclusively to syntaxial immobile vein systems; e.g., Bons et al., 2012), the viscosity of water at room temperature is = 10 −3 and much lower at higher temps. The final area (2D volume) is about 10 −6 − 10 −3 2 . The potential host lithologies are the same as for dikes, hence host rock fracture toughness and modulus are as above. Such a low viscosity fluid in a fracture of such small volume equilibrates almost instantly, and therefore it must be toughness-controlled. In large veins (>1 m aperture), complete sealing by mineral precipitation within the vein may occur very slowly or not at all (over years to millions of years; e.g., the calcite infills dated by Roberts & Walker, 2016) allowing ample time for the hydrofracture to relax toward equilibrium (if this was required). Hence veins are expected to be very strongly toughness-dominated ( → ∞ ).
In the context of kinetic versus toughness-controlled growth, it is also of interest to consider scaled analogue (laboratory) models, which each use different host materials and magma analogs. Here we focus on examples that aim to model dike ascent and materials that have measured values for fracture toughness: those that use a gelatin host analogue, typically with a low viscosity liquid (water or paraffin oil). For instance, using the constant flux experiments of , the injection rate (which is in 3 ∕ so converted here to 2D by dividing by the dike width, which is roughly the same as the height) is = 10 −5 2 ∕ and = 10 −4 − 10 −1 . The fracture toughness is not given, but can be determined from their Equation 7 and Figure 2 to be = 33 √ , and ′ = 10 3 . This gives 1 < < 10 . As with veins, this system is strongly toughness-dominated. Other analogue systems may fall outside of this range, such as those using granular mixtures (Schmeidel et al., 2017) or low-concentration laponite gels (Arachchige et al., 2021), and/or viscous fluids. The scaling mismatch in properties is noted elsewhere, in that for gelatine ′ = 10 −2 − 10 −1 (Kavanagh et al., 2013) whereas for rocks ′ = 10 −4 . This has the effect of increasing in those gelatine analogue models well into the fracture-controlled and equilibrium regime, and away from the region of natural dikes. Taking = 1 for veins in Figure 2 yields a prediction of ≈ 10 −3 √ for the host rock in these cases. Host rock data determined from laboratory tests for the dike and vein systems in Figure 1a suggests that ≈ 10 −4 √ in all cases. However, there is some evidence that modest increases in fracture toughness (by a factor of 3-5) could be possible at depth (Fialko & Rubin, 1997;Stoeckhert et al., 2016). Decreases in modulus are also possible at larger scales (Schultz, 1993) although this increase in compliance is largely due to the activity of joints which may be suppressed at depth. So we take = 10 −3 √ as the reference point. This is still consistent with the ranges assumed above, for example, = 1 √ and ′ = 1 . Plotting Equation 3 for values of ( ) on a thickness versus length diagram (Figure 3), veins exist at about ≈ 1 and dikes are about = 10 − 1, 000 . The exact results and position for dikes will therefore be dependent not only on the host rock properties, but also the magma flow rate and viscosity. Hence each dike system is unique and has the potential to occupy a different contour in ( ) . This finding becomes apparent on closer inspection of individual datasets for dikes in Figure 1a. Although the data are very scattered, power law fits are plotted with the = 0.5 exponent for all data (Olson, 2003). The Shiprock dikes ( = 0.44 ), Ethiopia dikes ( = 0.48 ) and Martian Elysium dikes ( = 0.43 ) each fit the = 0.5 model of Equation 3 reasonably well. On the other hand, the Karoo dikes give = 0.3 , the Sudan dikes = 0.22 , and the Deccan dikes = 0.06 , potentially indicating different conditions of emplacement for each set. In any case, dikes, veins, and analogue models, are not part of the same continuum and cannot be linked in these thickness versus length scaling plots.
Before progressing further, it is useful to reflect on the importance of the 2D nature of the model presented. Savitski and Detournay (2002) developed a higher order 3D model for kinetic-dominated growth in a penny-shaped crack increasing in volume over time as 3 = , where is a constant with units of 3 ∕ . Taking the length to be the dike diameter, the scaling relationship, in this case 7 of 13 gives a lower exponent of 0.25. However, this apparent conflict between the 2D and 3D models is easily resolved, as the final 2D model volume is 2 ≈ 4 and the final 3D model volume is 3 ≈ 4 3 .
2 8 meaning that = 2 3 , where is the final length of the dike. Substituting this scaling relationship in Equation 7 reproduces Equation 6 but with a different pre-factor (2.7) for the different geometry. The change in pre-factor is of little consequence, but this does raise the question about whether the magma injection rate depends on the final length of the dike, that is, by inference, the volume of magma emplaced. In part this would depend on whether the entire magma volume is available throughout dike growth, and/or how long can be physically sustained through magma supply. This is a question that cannot easily be answered, as it depends on many factors, such as the size of the magma packet that feeds the dike and its rate of ascent. However, it does not seem unreasonable to envisage that a large dike ( ≈ 100 ) might be fed somewhat more rapidly than a small dike ( ≈ 10 ) due to the enormous difference in the quantity of magma involved. The model of a size-invariant magma line source ( per meter) in this case appears to be more appropriate than a size-invariant point source ( ). As the former is more consistent with observations than the latter, we will proceed to develop the 2D model further, whilst noting that the 3D scaling can be obtained with the substitution of = 2 3 .
To summarize, the predictions of Figure 2 and the observations in Figures 1a and 3 both support the conclusion that dikes grow as non-equilibrium structures in the kinetic-dominated regime. Therefore, we now assume = 0 for the remainder of this paper. When considering non-equilibrium growth we propose that a dike extends in two phases: (1) an inflation phase, where the volume of magma in the dike increases over time; followed by (2) a relaxation phase, where the magma volume is fixed but the dike continues to extend, accommodated by magma flow, until it freezes. It is of interest to determine whether relaxation plays a significant role in dike scaling, but also to check that a dike cannot reach equilibrium within the predicted relaxation time. A similar model has been proposed previously for progression of horizontal sheet intrusions in − space, from (thick-short) laccolith to (thin-long) sill geometries (Bunger & Cruden, 2011) driven by magma body forces (the weight of the magma), but this does not apply in dikes.

Models for Non-Equilibrium Inflation and Relaxation Phases
In the context of the observed order of magnitude variation in the scaling relationship observations of Figure 1a, here we wish to develop a simple analytical solution which is a reasonable approximation of the full solution.
The work of Spence and Turcotte (1985) provides a good starting point for this. The novelty of the approach here is to extend their previous analysis for kinetics-dominated growth to allow a general expression for the volume evolution, 2 ( ) , such that non-linear inflation and relaxation can be considered, as illustrated in Figures 4 and 5. To model inflation, we assume power law growth with exponent , such that 2 ( ) = (Figure 4), then the − relationship of Equation 6 can now be written in a more general form (see Supporting Information S1) as where = 3 −1 3 +1 . Note that the exponent reduces to = 1 2 when = 1 and Equation 6 is recovered. To model the relaxation stage, we assume volumetric inflation ends at time = 0 with a final magma volume of 2 ( 0) = 0 = 0 . The dike can still evolve over time even without the addition of further magma, just at a much-reduced rate. If this evolution occurs for an additional relaxation time then the total time is = 0 + . In Supporting Information S1 we find that during inflation the length increases as 1 6 + 2 and during relaxation the length still increases but more slowly, tending toward 1 6 when ≫ 0 . Similarly, the thickness increases during inflation as − 1 6 + 2 and decreases during relaxation, tending toward − 1 6 when ≫ 0 . As such, lengthening due to relaxation occurs at its fastest immediately following inflation, and will slow rapidly (e.g., Figure 5). Even without accounting for the effect of cooling on viscosity, this means that a dike will not have sufficient time to reach equilibrium before it solidifies.
To estimate the relaxation time, we use the widely adopted model of Turcotte and Schubert (2002) for solidification of a dike: = 2 16 2 (9) Figure 4. Schematic illustration of linear dike growth (e.g., Spence & Turcotte, 1985) in the which the 2D volume (area) relates to the injection rate, , as a function of time, , relative to a non-linear magma injection model that uses a power law for the inflation stage up to time 0 , followed by an additional-constant volume-relaxation period of before the dike solidifies. Figure 5. An example inflation-relaxation sequence, showing the temporal evolution of dike length, , dike thickness, , dike pressure, , whereby a short, thick dike is rapidly injected over 4 days with an exponent of = 1 , leading to an increase in both its maximum (central) thickness and its length. At the end of the inflation phase, 0 , the dike relaxes the magma pressure over the following 16 days ( ) by a further increase in length, necessarily accommodated by a decrease in dike maximum thickness to conserve magma volume.

of 13
where is the Stefan constant and is the thermal diffusivity. Morita et al. (2006) calculated a value of = 0.36 and = 10 −6 2 −1 . This predicts a 1 m thick dike will take roughly 5.5 days to solidify, whereas a 5 m thick dike would take 140 days. We now have two different time scales: (a) inflation during magmatic volume increase ( 0 ); and (b) relaxation during constant magmatic volume ( ). As the dike thickness reduces during relaxation, such that ( ) , Equation 9 represent a quartic equation in terms of (see Equation C.2 in Supporting Information S1). Figure 6 illustrates a number of different dike growth and relaxation trajectories in − space. In Figure 6a we follow the observations of Morita et al. (2006) and take a magma injection rate = 1 2 ∕ , injection rate exponent = 0.65 (i.e., = 1 2 −0.65 ), plane strain modulus ′ = 10 , magma viscosity = 10 8 , and thermal diffusivity = 10 −6 2 ∕ . The blue line in Figure 6a, and subsequent plots, shows the inflation trajectory, with points along it showing the dike dimensions after different growth periods. For = 0.65 this has an exponent (slope) of = 0.33 . Once a dike stops increasing in volume, it progresses downward and to the right (increasing at the expense of ) along its relaxation trajectory (the dashed lines connecting points). This terminates in the green line, which signifies the end of the solidification time predicted by Equation 9. The green line represents an upper bound on the relaxation time, as it does not take into account cooling during growth, or any increase in viscosity during relaxation, though again it is noted that relaxation will be fastest when it starts. The full extent of relaxation is therefore hard to determine, but it is expected that a dike of a given volume will form somewhere between the blue line and the green line. In Figure 6a, it can be seen that these conditions envelop a significant portion of the observed dikes. A dike reaches the first point (a length of 59 m and a thickness of 2.8 m) in 20 min. The upper bound on its relaxation time is then about 4 days, substantially longer than the growth time. In this case, neglecting cooling during inflation is reasonable. As dikes get larger, the inflation time increases relative to the relaxation time, as the thickness is not increasing as rapidly as the volume. For the dike observed in Morita et al. (2006) inflation takes 8 days but relaxation could be as long as 100 days thereafter. The very largest dike shown in Figure 6a grows to a length of 63 km over 40 years, but with only a comparatively short 8 years to relax. In this case, neglecting cooling during inflation is not reasonable, but as the amount of relaxation undertaken is insignificant this is not important. The relaxation curve always has a higher exponent (slope) than the inflation line, and in this case the exponent increases to = 0.55 for the smallest dikes, converging back to = 0.33 for the largest dikes. Figure 6b shows that reducing the viscosity to = 10 6 drops the inflation and relaxation curves downwards, toward some of the thinner dike sets. In Figure 6c, a higher magma injection rate of = 10 2 −0.65 moves the growth curve upwards to encompass some of the thicker observed dikes, showing that the effect of relaxation could be quite substantial even in the larger length scales under this scenario. Figure 6d models rapid linear ( = 1 ) growth, for which nearly all the observed dikes sit between the inflation and relaxation curves. Figure 6e shows the effect of using a much lower exponent of = 0.5 (slower inflation), that is, insignificant relaxation for large dikes. Finally, Figure 6f shows that decreasing the thermal diffusivity to = 10 −7 2 ∕ leads to slower cooling and a wider zone between the inflation and relaxation curves as would be expected. Conversely, an increase in the thermal diffusivity will lead to a reduction in the zone of possibility for observed dikes.

Comparison With Natural Intrusions
In our model, dikes can extend their length in two stages: (a) an inflation stage in which both length and thickness increase, and (b) a constant volume relaxation stage, in which length can only grow at the expense of maximum thickness. In reality the relaxation stage is likely to be highly variable, and dependent on the details of the cooling and solidification processes. The model shown here assumes that cooling initiates at the onset of the second stage, whereas for major dikes it is much more likely that parts will cool during the initial stage of volume increase, due to contact conduction with the host rock walls. The temperature distribution within the magma will therefore be a minimum at the walls and increase to a maximum at the center. The picture is further complicated by the potential for temperature gain through the latent heat of crystallisation, and the increase in magma viscosity with crystal content. Here it is assumed that dike relaxation stops once it has solidified in the middle, at the position of maximum thickness, but this is not necessarily the case. Solidification within intrusions can be unevenly distributed, GILL ET AL.
10.1029/2022JB024120 10 of 13 leading to localisation of magma flow into channels (e.g., Holness & Humphreys, 2003). This localized flow of hot magma can lead to remobilization of accreted materials of variable viscosity across the conduit (e.g., Walker et al., 2017). Toward the tips, where the dike is much thinner, freezing could occur more rapidly. If the dike length is still extending at a sufficient rate (Delaney & Pollard, 1981) the magma at the tip will continue to be refreshed by an influx of hot material, preventing freezing. As such, the exact criteria that determines when a dike stops lengthening requires further investigation. The relaxation trajectories for length and thickness evolution shown as dashed green lines in Figure 6 therefore represent the maximum bound for relaxation. In nature then, some intermediary position of relaxation is probable, since a dike will undergo cooling during ascent, reheating as the magma crystallizes, and further cooling to an ambient geotherm, set within host rocks and accreted dike margins that have variable thermal diffusivity properties. Relaxation presents, therefore, an additional process that will result in − scatter for individual dikes within a larger volcanic system. The history of the freezing process will also determine the final internal pressure distribution, which will not be uniform or linear. This will be expressed in the final shape of the intrusion, which could result in a form between that of a lenticular geometry ; (c) higher growth rate = 10 2 −0.65 ; (d) higher growth exponent (and rate) with = 1 and = 1 (i.e., = 1 2 ∕ ); (e) lower growth exponent = 0.5 with increased = 10 2 −0.5 (note that without increasing , growth takes hundreds of years); (f) lower thermal diffusivity (which affects cooling rate) with = 10 −7 2 ∕ . Observed dikes are expected to lie in the region between the upper (solid blue) and lower (dashed green) lines under the stated conditions. The Shiprock dikes are shown separately here as they are individual echelon surface segments of a larger underlying dike (Scholz, 2010) and hence are not necessarily expected to comply with the model presented here. GILL ET AL.
10.1029/2022JB024120 11 of 13 with tapered tip profiles, and the elliptical to superelliptical profiles associated with an equilibrium pressure distribution shown in Figure 1b (Spence & Turcotte, 1985; see e.g., the schematic illustrations in Figure 5). This change in shape at the tip due to local magma redistribution in the final stages of freezing may change the stress distribution and failure mechanism at the tip Walker et al., 2021), which may affect intrusion lengthening, thereby introducing further scatter in − space.
The question remains whether relaxation is evident during active intrusion and within the rock record. In active systems, Morita et al. (2006) provide some evidence for two stage growth in their study of earthquake swarms during dike intrusion in Izu Peninsula, Japan. From geodetic observations, they show that volume increase during dike growth occurred over 14 days, whereas associated seismicity occurred for 20 days. Based on their dimensions, the relaxation model here would have a conservative prediction for relaxation on the order of 100 days (maximum 400 days), which is far in excess of the 6 days indicated in the Morita et al. (2006) study. There are two immediate explanations for the discrepancy: (a) our model is an overestimate because there is likely to be significant cooling and potential solidification during the volumetric growth of the dike; and (b) fracture growth to accommodate relaxation may fall below seismicity detection limits (i.e., it becomes aseismic), particularly if growth is accommodated by dominantly tensile failure in the host rock (i.e., non-double couple mechanisms) as opposed to the shear-fracturing (double-couple mechanisms) shown in most dike seismicity studies. The latter explanation appears to be the case even for the volumetric inflation stages elsewhere, such as the Bárðabunga-Holuhran diking event in Iceland (Sigmundsson et al., 2015). Emplacement of the Bárðabunga-Holuhran dike induced earthquakes during growth laterally and toward the surface for about 2 weeks (Ágústsdóttir et al., 2016), followed by a 6-month eruption phase, and a further 6 months of post eruption seismicity along the length of the dike section (Woods et al., 2019). Ágústsdóttir et al. (2016) interpret post-eruption earthquakes detected at 5-7 km depth as representing late-stage equilibration of magma pressure in the dike; that is, relaxation. Geodetic measurements indicate that seismicity did not capture all pre-eruptive dike growth at shallow depths, including that necessary for magma to reach the surface (Sigmundsson et al., 2015), hence it is conceivable that some late stage growth may also go undetected in the shallow crust, particularly where dike lengthening is accommodated by tensile failure of the host rock (Ágústsdóttir et al., 2016;Rubin et al., 1998) or by dilatation of existing structures . In any case, it is worth noting that eruption will have served to reduce excess magma pressure in the remaining dike, in which case the actual period of relaxation should be greatly reduced compared to the timescales predicted in Figure 6. In addition, the formation of a graben above the dike would likely place further constraints on the dike's ability to relax, as thickness reduction could require reactivation, and potentially inversion, of the graben fault system.
Lengthening during the relaxation stage may be a cryptic feature in the rock record also, since the diagnostic feature of relaxation is lengthening without volume increase, requiring thinning at some position along the dike (e.g., Daniels et al., 2012). Field-based studies of frozen intrusions have shown the potential for late-stage lengthening at preserved tip zones, that can be identified from overprinting textures and tip zone deformations Walker et al., 2021). The tip forms of such intrusions are typically blunted, with squared-ends and a relatively constant thickness compared to the bladed geometry that should result from rock splitting. In a linear elastic framework, this constant thickness would represent a constant magma pressure. However, such examples are commonly associated with distributed shear faulting at the tip, within the intruded host rock, which is interpreted to represent the magma front moving forward as a viscous indentor (Galland et al., 2019;Spacapan et al., 2017). This may still represent a constant pressure in the conduit, but could represent a plug flow of relatively cool and high viscosity magma. As noted above, introducing a plug flow regime is equivalent to changing the effective channel thickness in the model, and would therefore influence the relaxation model. Viscous inden tation is also a relatively inefficient growth mechanism, particularly compared to elastic (tensile) splitting of the host rock, since the newly created fracture surfaces remain in contact and maintain a residual friction. Although a dike may grow by this mechanism over short distances, it is also possible that residual magma pressure may activate a new and more efficient pathway elsewhere on the dike , leading to only very local lengthening, and reducing the likelihood of observing the true maximum length dimension of the dike. In any case, these features are not necessarily uniquely related to a relaxation stage of growth, and further study would be required to constrain the distribution of such features at the periphery of individual dikes relative to changes in the thickness. 12 of 13

Conclusions
Toughness-dominated models for dike growth predict unreasonably large values for the rock fracture toughness, based on the assumption that magma pressure is constant within the dike, despite the need for pressure gradients to drive magma flow. However, this approach has attracted most attention in the literature as it describes a final equilibrium (static) state. We address this problem using a kinetic-dominated analytical approach to consider 2D dike growth and geometry. Dike growth in a kinetic analysis can be split into two stages, with an inflation stage characterized by (horizontal) lengthening and thickening as the magma volume increases, followed by a relaxation stage in which the magma volume is fixed and pressure gradients are relieved within the dike. Conservation of volume during this second stage of growth requires that increases in dike length must be met with a decrease in maximum thickness. Lengthening is shown to exhibit a 1 6 scaling over time, hence it slows quite rapidly, allowing ample time for complete cessation of dike growth due to solidification of the magma. We find that altering the conditions during dike growth (such as the magma injection rate and viscosity, and the host rock plane strain modulus and thermal diffusivity) within a reasonable range possible in nature, results in a range of potential pathways to a final thickness to length scaling ratio. Although observed exhumed dikes generally show a = 0.5 relationship (ultimately with significant scatter), there is no unique route by which they can achieve their final scaling.

Data Availability Statement
Data presented in this paper is available via the University of Leicester Figshare portal at https://doi.org/10.25392/ leicester.data.20559855.