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.

Please refer to the supplement