Models of buoyancy-driven dykes using continuum plasticity and fracture mechanics: a comparison

Li, Yuan; Davis, Timothy; Pusok, Adina E.; Katz, Richard F.

doi:https://doi.org/10.5194/egusphere-2024-3504

Preprints

https://doi.org/10.5194/egusphere-2024-3504

Preprints

05 Dec 2024

| 05 Dec 2024

Models of buoyancy-driven dykes using continuum plasticity and fracture mechanics: a comparison

Yuan Li, Timothy Davis, Adina E. Pusok, and Richard F. Katz

Abstract. Magmatic dykes are thought to play an important role in the thermomechanics of tectonic rifting of the lithosphere. Our understanding of this role is limited by the lack of models that consistently capture the interaction between magmatism, including dyking, and tectonic deformation. While linear elastic fracture mechanics (LEFM) has provided a basis for understanding the mechanics of dykes, it is difficult to consistently incorporate LEFM into geodynamic models. Here we further develop a continuum theory that represents dykes as plastic tensile failure in a two-phase, Stokes–Darcy model with a poro- viscoelastic–viscoplastic (poro-VEVP) rheological law (Li et al., 2023). We validate this approach by making quantitative comparison with LEFM, enabled by a novel poro-LEFM formulation. The comparison shows that dykes in our continuum theory propagate slowly—a consequence of Darcian drag on the magma. Moreover, dissipation of mechanical energy in the poro-VEVP model implies a high critical stress intensity in LEFM. We improve the poro-VEVP model by reformulating the compaction stress and incorporating anisotropic permeability in regions of plastic failure.

Received: 09 Nov 2024 – Discussion started: 05 Dec 2024

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this preprint. The responsibility to include appropriate place names lies with the authors.

Download & links

Yuan Li, Timothy Davis, Adina E. Pusok, and Richard F. Katz

Status: final response (author comments only)

CC1:
'Comment on egusphere-2024-3504', Giacomo Medici, 11 Dec 2024

General comments
Very good multidisciplinary research with a focus on dykes, I definitely enjoyed reading it. The amount of literature on the hydraulic properties of fractured rocks at such a large scale is not large. Please, see my specific comments to improve the manuscript.

Specific comments
Lines 62-70. Porosity. Total or effective porosity? When I think to the porosity of fractured rocks, I tend to consider this concept. Think if you need to specify something
Lines 10, 71-76. Permeability anisotropy. Vertical or horizontal? You need to explain this point in the abstract.
Line 73. Permeability tensor. Are you thinking to apply this concept at which observation scale?
Line 73. Permeability tensor. If you think that is a good idea to discuss the observation scale, please provide details on the depth and lateral extension.
Lines 74-75 “Anisotropic permeability can arise from anisotropic stresses and aligned pores or fractures”. Please, insert recent and relevant literature on the anisotropic permeabilities due to either anisotropic stress and orientation of fractures:
- Medici G, Ling F, Shang J 2023. Review of discrete fracture network characterization for geothermal energy extraction. Frontiers in Earth Science, 11, 1328397.
- Lei, Q., Latham, J.P. and Tsang, C.F., 2017. The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks. Computers and Geotechnics, 85, 151-176.

Lines 174. You reference multiple times Snow 1979. Do you need to discuss the hydraulic / mechanical aperture of the joints. What about the cubic law?

Figures and tables
Figure 1. Do you need to insert a spatial scale in your conceptual model?
Figure 3a. The figure describes the porosity and solid deformation field at t = 2 kyr. This is an important figure and if the reader wants to catch the details need to zoom in a lot. Please, enlarge the size
Figure 3a. Can this figure be separated from the others?
Figure F1. This is a very important figure from a conceptual point of view. Indeed, the image shows physical variations as a function of the depth. If you introduce this figure in the main body
of the manuscript, you would rise either the readability or the impact of your research.

Citation: https://doi.org/10.5194/egusphere-2024-3504-CC1
- AC1: 'Reply on CC1', Yuan Li, 13 Jan 2025
  
  We thank the reviewer for this feedback. Below, we address each specific comment in turn, with our responses indicated by R:.
  * Lines 62-70. Porosity. Total or effective porosity? When I think to the porosity of fractured rocks, I tend to consider this concept. Think if you need to specify something.
  R: In this manuscript, we assume total and effective porosity are equivalent, implying that all voids within the grain matrix are interconnected and fully saturated with liquid. We acknowledge that a distinction between these porosities exists and is important in some contexts. However, for partially molten rock in the asthenosphere, the pore structure is controlled by textural equilibrium, which results in a highly connected pore network. Therefore, equating the total and effective porosity is common practice in the geodynamic modelling community for the scales and processes considered. Thus, further elaboration is not included here.
  * Lines 10, 71-76. Permeability anisotropy. Vertical or horizontal? You need to explain this point in the abstract.
  R: To clarify, we will revise the sentence in line 73 as follows:
  Original: ``We resolve this discrepancy by introducing an anisotropic permeability tensor into the poro-VEVP model...”
  Revised: ``We resolve this discrepancy by introducing an anisotropic permeability tensor into the two-dimensional poro-VEVP model...”
  This revision makes clear that the anisotropic permeability is implemented within a 2D model, which is later described as occupying a vertical plane.
  * Line 73. Permeability tensor. Are you thinking to apply this concept at which observation scale?
  R: The ``observation scale" is not a defined quantity in the context of our mathematical model. We define anisotropic permeability at the subgrid scale, where it is described by the tensor. Hence the anisotropy is uniform within a grid cell; anisotropy at larger scales emerges as part of the numerical solution. This large scale structure might be compared with observations.
  * Line 73. Permeability tensor. If you think that is a good idea to discuss the observation scale, please provide details on the depth and lateral extension.
  R: Thank you for the suggestion. However, we opt not to discuss this topic, as it is less relevant to the manuscript’s main focus.
  * Lines 74-75 “Anisotropic permeability can arise from anisotropic stresses and aligned pores or fractures”. Please, insert recent and relevant literature on the anisotropic permeabilities due to either anisotropic stress and orientation of fractures:
  - Medici G, Ling F, Shang J 2023. Review of discrete fracture network characterization for geothermal energy extraction. Frontiers in Earth Science, 11, 1328397.
  - Lei, Q., Latham, J.P. and Tsang, C.F., 2017. The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks. Computers and Geotechnics, 85, 151-176.
  R: Thank you for suggesting these relevant studies. We will include them in the references of the revised manuscript.
  * Lines 174. You reference multiple times Snow 1979. Do you need to discuss the hydraulic / mechanical aperture of the joints. What about the cubic law?
  R: We do not need to discuss the hydraulic aperture of the joints and its effect on permeability models as Snow (1979) did. This is a subgrid feature that we cannot resolve. Such a detailed discussion lies beyond the scope of this manuscript, which focuses on comparing the continuum (poro-VEVP) and discrete (LEFM) models given the same constitutive laws.
  * Figure 1. Do you need to insert a spatial scale in your conceptual model?
  R: We appreciate the suggestion but believe a spatial scale is unnecessary for this figure. The schematic is designed to highlight differences between the classical LEFM model and the proposed poro-LEFM model. For such a conceptual diagram, it is not necessary to either show dimensional sizes or the relative size ratio, such as the head size to the far-field aperture. A similar sketch is commonly used in LEFM studies, for example Figure 1 in Roper and Lister (2007).
  - Roper, S. M. and Lister, J. R., 2007. Buoyancy-Driven Crack Propagation: The Limit of Large Fracture Toughness, Journal of Fluid Mechanics, 580, 359–380.
  * Figure 3a. The figure describes the porosity and solid deformation field at t = 2 kyr. This is an important figure and if the reader wants to catch the details need to zoom in a lot. Please, enlarge the size.
  Figure 3a. Can this figure be separated from the others?
  R: Thank you for this suggestion to improve the presentation of results. We will separate Figure 3a from other panels to make it bigger.
  * Figure F1. This is a very important figure from a conceptual point of view. Indeed, the image shows physical variations as a function of the depth. If you introduce this figure in the main body of the manuscript, you would rise either the readability or the impact of your research.
  R: We appreciate this suggestion. However, after careful consideration, we believe it is more appropriate to leave Figure F1 in the Appendix. This figure highlights a complex difference that, while interesting, has limited impact on the main conclusions and is not yet fully understood. As such, we include it as supplementary material for specialist readers without diverting focus from the primary objectives of the manuscript.
  
  Citation: https://doi.org/10.5194/egusphere-2024-3504-AC1
RC1:
'Comment on egusphere-2024-3504', Anonymous Referee #1, 27 Jan 2025

The manuscript by Li et al. discusses the possibility of modeling dykes in their in-house code for the two-phase Stokes–Darcy model with a poroviscoelastic–viscoplastic rheological law. Li et al. implemented dyke with the hyperbolic yield criterion assuming strong anisotropy of permeability inside a dyke. The authors further compare their results with the LEFM model by modifying the latter by introducing porous Darcy flow inside a dyke. While the authors show the vertical high-porosity region produced by their model, I wonder if their results are relevant to dyke propagation processes. First, they had to introduce a strong anisotropy of permeability inside a channel they produced, and the flow was essentially driven by buoyancy and this anisotropy. Thus, dyke always follows the direction of the anisotropy. However, most of the rocks have higher permeability in the horizontal direction, and the planes of weakness of anisotropic rocks are also in the horizontal direction. Thus, this argument would be more suitable to simulate sills. In reality, sills and dykes often accompany each other forming interconnected plumbing systems. Second, they “prescribe” dyke by setting the initial elevated porosity region at the bottom of their computational domain prescribed by anisotropic Gaussian porosity distribution. And third, the width of their dyke is one grid size. This would suggest that at higher resolution, their dyke will be an infinitely thin line. In reality, dykes always have a finite length, sometimes meters. The authors seem unaware of the state-of-the-art observations and models. The introduction contains a review of their own work, ignoring the work of other groups. In addition, the model equations presented in the paper suffer from a lack of mathematical rigor. The mathematical errors must be fixed to assess the validity of their results. Thus, I am not convinced by the described model. See detailed comments below:

Line 12. “Magmatic dykes, formed by fluid-driven fracture, are important pathways for magma ascent across the lithosphere.” This is just one of the hypotheses; correct the sentence and cite other possible mechanisms.
Line 22. “In most previous work, the mechanics of dykes is formulated in terms of linear elastic fracture mechanics (LEFM).” This is not entirely true. The models involving plastic failure as a mechanism for dyking were presented in (e.g., Alkhimenkov et al., 2024; Gerbault et al., 2018; Gerbault et al., 1998; Minakov et al., 2017).
Line 24. “LEFM conceptualises dykes as mode-I fractures opened at the tip and widened by magma flow (Rivalta et al., 2015)” There are earlier than that references suggesting mode-I fracture as a dyking mechanism (Mckenzie et al., 1992; Ode, 1957).
Lines 35-49. This is a very biased review. The statement “In the geodynamic context of the hot, ductile asthenosphere, magma transport has long been modelled using a poro-viscous, Stokes–Darcy theory … These studies were limited to hot asthenospheric regions by the use of a purely viscous rheological law” is wrong. Viscoelastic and viscoelastoplastic models were also used in geodynamics (e.g., Bercovici et al., 2001; Burov, 2011; Cai and Bercovici, 2013; Connolly and Podladchikov, 1998; Kaus and Podladchikov, 2006; Kiss et al., 2023; Petrini et al., 2020).
Lines 47-48. “… and (de)compaction in poro-viscous dynamics (Katz et al., 2022)…” Strange reference to “(de)compaction”. Include a reference to the original paper introducing this term (Yarushina and Podladchikov, 2015).
Lines 57-70. This section lacks a review of the models of elastic fracture with cohesive zones. Such a model was first proposed by Barenblatt (1962). The process zone concept at the fracture tip has been actively used ever since (e.g., Viesca et al., 2008). How does your estimation of fracture toughness and energy compare with the estimation of these parameters in the classical Barenblatt’s model?
Lines 71-74. “Isotropic permeability within the poro-VEVP dyke promotes widening by horizontal porous flow, a behaviour not associated with real (or LEFM) dykes. We resolve this discrepancy by introducing an anisotropic permeability tensor into the poro-VEVP model to limit leakage and enhance fracture propagation…” Discuss other possible reasons for not reproducing dyke. For example, an underpressure at the fracture tip generated at the moment of fracture opening can be a driving force for the flow. Fracture opening must be must faster process than Darcy flow, and thus, there would be a region with undrained conditions at the tip where pore volume is increased due to fracture generation that would lead to the pressure drop at the tip, increasing pressure gradient and thus the flow.
Line 105. This equation is not standard, give explanations.
Line 108. Why do you assume fluid pressure is equal to solid pressure? Further in the text you assume that p_l = const and arrive at equation (5), while p_s given by equations (2) and (3) is not constant. These two assumptions are not consistent.
Line 147. Equation (6) and the corresponding equation in the Appendix have a wrong gravity term. The conservation of total momentum has a total density in the gravity term (Bercovici et al., 2001; Yarushina and Podladchikov, 2015).
Line 148. “Here, (1−ϕ)t^s and −(1−ϕ)ΔP represent the effective shear and decompaction stresses” The term (1−ϕ)ΔP=P_total - P_l = P_e is commonly referred to as Terzaghi’s effective pressure (Coussy, 2004; Detournay and Cheng, 1993; Nur and Byerlee, 1971; Rice and Cleary, 1976; von Terzaghi, 1923). There is no need to invent new names, which might confuse the reader. Furthermore, not recognizing common concepts might lead to further mistakes. There is a further lack of consistency in the notations in the appendix.
Lines 156 – 159. This is already an algorithmic formulation. Please provide the governing equations you are solving. These are absent in the main text and in the appendix.
Line 160. “Here, C is the solid decompaction rate…” Again, you are introducing non-common terms for very standard parameters. Looking at the appendix, I might assume that C is the viscous or inelastic part of the volumetric strain rate. Using conventional terminology would help in understanding your model to a broader audience and avoiding mistakes.
Lines 169-170. “This approach only changes how the stress-balance equation is linked with the plastic yield condition, without altering either of the two physics.” Explanation is confusing at best. I do not see how physics stays the same.
Line 209. Since you do not provide your rheology equations where G and Z are parameters, it is useless to provide these expressions here.
Line 214. “Appendix C reviews the full system of equations.” This is not true. This appendix contains algorithmic equations implemented in the code that include purely numerical parameters together with physical quantities, not the full system of governing equations. Still, a similar statement should appear already in line 149.
Line 220. It would be helpful with the 2D figure illustrating the initial distribution of porosity.
Line 248. “Appendix D provides details of the formulation for each local work rate.” Incorrect. You provide this estimation only for one in the appendix.
Line 255. “Then, assuming that the elastic contributions to these work rates are approximately equal,” What is the basis for this assumption?
Line 260. Why propagation speed is assumed constant?
Line 292. “the tip advances approximately the same distance in each 0.4 kyr interval” Why? Discuss mechanisms and reasons for such behavior.
Line 296-299. The 2D figure of plastic failure around the tip would help in understanding. You are referring everywhere in the text that your dyke is generated by mode I fracture. How can you be sure that it is mode I and not mode II mode, given that your yield criterion contains both? This statement needs further proof and illustration. Are there any shear stresses associated with dyke propagation?
Lines 302-303. “In the poro-VEVP model, buoyancy induces plastic tensile failure throughout the head region, whereas in the poro-LEFM model, fracture…” Comparison with the Barenblatt model, not just LEFM, would be appropriate here.
Figure 5 caption. “the liquid overpressure ΔP.” Yet, one more term is invented here for the same parameter (see my comment on line 148). Consider using conventional effective pressure.
Lines 372-373. “The dyke width is determined by the grid size, which is a limitation of the present discretized solutions of the continuum models.” Discuss these limitations and ways to overcome them. The problem of mesh dependency and respective regularization is widely discussed in the geodynamic community (e.g., Duretz et al., 2023).
Lines 418-422. Maybe anisotropy is not the right mechanism? The results do not look convincing.
Lines 470-480. This description does not add clarity. Consider removing it and adding appropriate references in the main text.
Line 525. These are algorithmic equations; what are the physical relations?
Line 526. What is overpressure and dilatancy pressure here? Consider using conventional terminology. You already used equation B5 when deriving C1 from (6). You need a separate independent equation for it.
Lines 530-532. Again, comment on terminology. You introduce here effective pressure, which had to be done much earlier. However, your definition of P_e is inconsistent with commonly used definitions for effective pressure in elastoplastic materials (see comment on line 148). The same P_e goes into what seems to be the “standard” elastoplastic flow rule given by equation (C7). Such a choice of P_e might cause problems with the thermodynamic admissibility of your system of equations. Coussy (2004) shows that for materials with plastically incompressible solid grains, the effective pressure used in the flow rule should be Pe = Ptotal -Pl. You also assumed the incompressibility of grains (line 111). Thus, I have serious doubts about the mathematical correctness and eligibility of the presented model.
Lines 539-541. “For generality…” Equations (C9) are not general. Where do they come from?
Lines 546 – 580. This whole analysis is meaningless here as the full set of equations is not given.
Lines 551-552. “The local work rate associated with deformation at a point can be expressed as the product of the strain rates and effective stresses causing the deformation (Batchelor, 2000; Katz, 2022).” The concept of work rate is much more fundamental and dates much earlier than 2000 or 2022. Consider more appropriate references here.
Line 555. Missing 1/3 multiplier in the equation. C seems to be a viscous part, not the full strain rate. But again, the missing governing equations do not allow final judgment.
Line 570-576. Equations (D7) and (D8) need further explanations. You assumed that for viscoplastic stresses. However, this is not what is written as a viscoplastic constitutive equation (C7). Provide the proper derivation.

References:
Alkhimenkov, Y., Khakimova, L., Podladchikov, Y., 2024. Shear Bands Triggered by Solitary Porosity Waves in Deforming Fluid-Saturated Porous Media. Geophysical Research Letters 51.
Barenblatt, G.I., 1962. The Mathematical Theory of Equilibrium Cracks in Brittle Fracture, in: H.L. Dryden, T.v.K., G. Kuerti, F.H. van den Dungen, L. Howarth (Ed.), Advances in Applied Mechanics. Elsevier, pp. 55-129.
Bercovici, D., Ricard, Y., Schubert, G., 2001. A two-phase model for compaction and damage 1. General Theory. J Geophys Res-Sol Ea 106, 8887-8906.
Burov, E.B., 2011. Rheology and strength of the lithosphere. Mar Petrol Geol 28, 1402-1443.
Cai, Z.Y., Bercovici, D., 2013. Two-phase damage models of magma-fracturing. Earth Planet Sc Lett 368, 1-8.
Connolly, J.A.D., Podladchikov, Y.Y., 1998. Compaction-driven fluid flow in viscoelastic rock. Geodin Acta 11, 55-84.
Coussy, O., 2004. Poromechanics, 2nd ed. Wiley, Chichester, England; Hoboken, NJ.
Detournay, E., Cheng, A.H.D., 1993. Fundamentals of poroelasticity, in: Hudson, J.A. (Ed.), Comprehensive rock engineering : principles, practice and projects. Pergamon Press, Oxford, pp. 113-171.
Duretz, T., Räess, L., de Borst, R., Hageman, T., 2023. A Comparison of Plasticity Regularization Approaches for Geodynamic Modeling. Geochem Geophy Geosy 24.
Gerbault, M., Hassani, R., Lizama, C.N., Souche, A., 2018. Three-Dimensional Failure Patterns Around an Inflating Magmatic Chamber. Geochem Geophy Geosy 19, 749-771.
Gerbault, M., Poliakov, A.N.B., Daignieres, M., 1998. Prediction of faulting from the theories of elasticity and plasticity: what are the limits? J Struct Geol 20, 301-320.
Kaus, B.J.P., Podladchikov, Y.Y., 2006. Initiation of localized shear zones in viscoelastoplastic rocks. J Geophys Res-Sol Ea 111.
Kiss, D., Moulas, E., Kaus, B.J.P., Spang, A., 2023. Decompression and Fracturing Caused by Magmatically Induced Thermal Stresses. J Geophys Res-Sol Ea 128.
Mckenzie, D., Mckenzie, J.M., Saunders, R.S., 1992. Dike Emplacement on Venus and on Earth. J Geophys Res-Planet 97, 15977-15990.
Minakov, A.V., Yarushina, V.M., Faleide, J.I., Krupnova, N., Sakoulina, T., Dergunov, N., Glebovsky, V., 2017. Dyke emplacement and crustal structure within a continental large igneous province - northern Barents Sea. Geological Society, London, Special Publications.
Nur, A., Byerlee, J.D., 1971. Exact effective stress law for elastic deformation of rock with fluids. Journal of Geophysical Research 76, 6414-6419.
Ode, H., 1957. Mechanical Analysis of the Dike Pattern of the Spanish Peaks Area, Colorado. Geol Soc Am Bull 68, 567-575.
Petrini, C., Gerya, T.V., Yarushina, V.M., van Dinther, Y., Connolly, J.A.D., Madonna, C., 2020. Seismo-hydro-mechanical modelling of the seismic cycle: methodology and implications for subduction zone seismicity. Tectonophysics, 228504.
Rice, J.R., Cleary, M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Reviews of Geophysics 14, 227-241.
Viesca, R.C., Templeton, E.L., Rice, J.R., 2008. Off-fault plasticity and earthquake rupture dynamics: 2. Effects of fluid saturation. J Geophys Res-Sol Ea 113.
von Terzaghi, K., 1923. Die berechnung der durchlassigkeit des tones aus dem verlauf der hydromechanischen spannungserscheinungen. Sitzungsbericht der Akademie der Wissenschaften (Wien): Mathematisch-Naturwissenschaftlichen Klasse 132, 125-138.
Yarushina, V.M., Podladchikov, Y.Y., 2015. (De)compaction of porous viscoelastoplastic media: Model formulation. Journal of Geophysical Research Solid Earth 120.

Citation: https://doi.org/10.5194/egusphere-2024-3504-RC1
- AC2: 'Reply on RC1', Yuan Li, 14 Mar 2025
  
  The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-3504/egusphere-2024-3504-AC2-supplement.pdf
  
  Citation: https://doi.org/10.5194/egusphere-2024-3504-AC2
RC2:
'Comment on egusphere-2024-3504', Anonymous Referee #2, 07 Feb 2025

Please refer to the supplement

Citation: https://doi.org/10.5194/egusphere-2024-3504-RC2
- AC3: 'Reply on RC2', Yuan Li, 14 Mar 2025
  
  The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-3504/egusphere-2024-3504-AC3-supplement.pdf
  
  Citation: https://doi.org/10.5194/egusphere-2024-3504-AC3

Yuan Li, Timothy Davis, Adina E. Pusok, and Richard F. Katz

Viewed

Total article views: 292 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
210	69	13	292	13	11

HTML: 210
PDF: 69
XML: 13
Total: 292
BibTeX: 13
EndNote: 11

Views and downloads (calculated since 05 Dec 2024)

Month	HTML	PDF	XML	Total
Dec 2024	91	23	5	119
Jan 2025	54	20	4	78
Feb 2025	41	10	1	52
Mar 2025	15	7	0	22
Apr 2025	9	9	3	21

Cumulative views and downloads (calculated since 05 Dec 2024)

Month	HTML	PDF	XML	Total
Dec 2024	91	23	5	119
Jan 2025	54	20	4	78
Feb 2025	41	10	1	52
Mar 2025	15	7	0	22
Apr 2025	9	9	3	21

Viewed (geographical distribution)

Total article views: 298 (including HTML, PDF, and XML) Thereof 298 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 23 Apr 2025

Short summary

Magmatic dykes transport magma to the Earth's surface, sometimes causing eruptions. We advanced a model of dyking, treating it as plastic deformation in a porous medium, unlike the classic model that treats dykes as fractures in elastic solids. Comparing the two, we found the plastic model aligns with the fracture model in dyke speed and energy consumption, despite quantitative differences. This new method could be a powerful tool for understanding volcanic processes during tectonic activity.


Total:	0
HTML:	0
PDF:	0
XML:	0