General comments

This is a well written and organised paper that considers the effects of increased permeability in shallow-level igneous sills on hydrothermal fluid flow and both the formation and transport of hydrocarbons. The study builds on field observations of cooling joints infilled with bitumen, as well as dm-thick bitumen dykes in the surrounding host rock of igneous sills in the Neuquén Basin, Argentina. Raman spectroscopy of two bitumen samples collected from veins at the outer intrusion margin indicates hydrocarbon temperatures of 350-500 ºC while circulating through the sill. These findings in combination with field observations of fractured intrusions were used to design a series of conceptual FEM/FVM fluid flow simulations to explore how increasing permeability within the cooling intrusion caused by the formation of cooling joints affects hydrothermal fluid flow and the transport of hydrocarbons. In their modelling results, the authors identify three distinct flow phases of fluids and hydrocarbons around a sill, which they use to discuss the evolution of the igneous petroleum systems in the northern Neuquén Basin. The main contribution of this study is the implementation of a simplified dynamic permeability to mimic cooling joints, which has a clear impact on hydrothermal fluid flow patterns. These findings are not only applicable to the transport of hydrocarbons, but also to sills in hydrothermal systems in general. Therefore, this study is a valuable contribution to the ongoing investigation of fluid flow in magmatic hydrothermal systems. However, I have concerns that the permeabilities used in the simulations are not representative of the study area which makes me wonder if the presented results can be used to interpret implications for igneous petroleum systems in the Río Grande Valley. Please see more detailed comments below. I would further like to request the authors to consider a couple clarification to: (1) their description of the used model and the model set-up, and (2) their discussions.

Specific comments

1) Permeability model of the host rock

I am not convinced that the described porosity-dependent permeability model for the host rock is reasonable to represent a low-permeability shale. Shales are commonly described as a low-permeability rock; e.g., Goral et al., 2020 (doi.org/10.1038/s41598-019-56885-y) state a maximum permeability typically less than 1000 nD (~9.8E-19 m²), which would not allow for significant advection (e.g., Ingebritsen et al., 2010; doi.org/10.1029/2009RG000287).

If there is no field evidence for these highly permeable shales in the Neuquén Basin, the presented models are not appropriate to interpret the igneous petroleum system and hydrocarbon transport in the Río Grande Valley. In that case, I would like to request the authors to either (1) refer to the host rock as a more permeable rock type (e.g., sandstone) that matches a permeability of ~1E-16 – 1E-14 m², or (2) to test/proof if simulations with host rock permeabilities of < ~1E-18 m² would give similar results as presented and discussed in this contribution. In the prior case, the changed host rock type should be considered in the discussion. E.g., What are the effects on organic matter transformation to methane compared to shale? Due to the change in host rock, the presented flow simulations would further not represent the field locality such that an interpretation of the implications for igneous petroleum systems in the Río Grande Valley may not be feasible.

While reading the section on host rock permeability, I was wondering why no brittle-ductile transition (BDT) was implemented for the host rock? Other studies on hydrothermal fluid flow suggest a BDT starting at temperatures of 360 ºC (e.g., Hayba and Ingebritsen, 1997; doi.org/10.1029/97JB00552). Would a BDT in the host rock decrease the permeability within the high-temperature aureole and thus affect fluid flow pathways?

2) Permeability model of the intrusion

Please provide more information on how the dynamic permeability within the intrusion is calculated. At the moment it is a bit unclear to me.

For your setup, I understand that the intrusion is impermeable (1E-20 m²) at T>=1100 ºC. Using the described linearised, temperature-dependent definition of the melt fraction, a crystallinity of 50% is reached at T=1000 ºC, which is the BDT and defines the onset of fracturing due to cooling. The next step is not clear to me. Does the permeability linearly increase during cooling until T=900 ºC is reached where the intrusion reaches the maximum permeability (1E-15 m²)? Is that correct? If so, the authors could refer to Iyer et al. (2013) who used a similar linear permeability approximation in one of their models.

Including the permeabilities of the intrusion for the threshold temperatures (T=1000 ºC, T=900 ºC) within the manuscript could also help to clarify the permeability model.

I would further like to invite the authors to justify and discuss the permeability values chosen for the fractured intrusion (1E-15 m²). Measured permeabilities of fractured intrusions within the Neuquén Basin (Spacapan et al., 2020) indicate permeabilities of ~5E-18 to 5E-15 m², with the majority of the samples being <5E-16 m². In their models, the authors use a permeability of 1E-15 m² for a fully solidified and fractured intrusion. Although this is only slightly above the maximum permeability reported by Spacapan et al. (2020), I would like to invite the authors to discuss potential effects of lower permeabilities as observed within intrusions in the Neuquén Basin on fluid flow and hydrocarbon transport. Would there be a fluid flow Phase 2 (“flushing”) also for permeabilities of <=1E-16 m² ?

As discussed in previous studies, permeability is a key controlling parameter within hydrothermal systems; e.g., the limiting permeability that allows for significant heat advection is 1E-16 m² (e.g., Ingebritsen et al., 2010). Therefore, it is critical to carefully decide on (and justify) the permeabilities used for both the host rock and the solidified intrusion as they will control the dynamics of the whole hydrothermal system including maturation and methane transport.

3) Model description

The authors refer to Galerne and Hasenclever (2019) when describing the model. However, a more detailed description of the model including the governing equations would be beneficial, given that this is the focus of the manuscript. This section could also be provided as appendix.

In addition, simplifications and assumptions of the model setup could be introduced and justified here (Section 3.1). The model considers a single-phase flow of a compressible fluid following Darcy’s law (L 251-252). From Galerne and Hasenclever (2019), I learned: “Throughout the calculations, pore pressures are above the critical point of pure water so that the fluid remains in a single‐phase state. Our single‐phase hydrothermal model requires this assumption, because in the system H2O–NaCl–CO2–CH4 phase transitions would be possible even at higher pressures.”. This assumption should be included and explained in the presented manuscript, and potential effects on the modelling results should be considered and discussed later in the manuscript. Given the shallow emplacement depth of 1-3 km, pore fluid boiling is plausible. Would this phase transition change fluid flow patterns and the transport and accumulation of the hydrocarbons?

Buoyancy effects caused by methane in the fluid are not considered in the models (L 262-264). Is this because the effects are too minor to significantly affect fluid flow or would buoyancy change the modelling results?

The exact domain size of the models is not given and should be included in the model description. It would also be important to explain/show how far away from the sill tip the no flow boundary is located.

Is the element size of 0.5 m consistent throughout the whole domain or only in a refined area within and around the intrusion?

4) Results and Figures

The figures are of good quality and visualise the results described in the manuscript. However, not all data discussed in the text are presented in the figures (Fig. 6, 7, 8, 9). Fluid pressure and the pore pressure distribution are both used to explain changes within the distinct fluid flow phases (i.e., “flushing” and “post-flushing”) and are also used throughout the discussion to explain fluid flow and methane transport. Therefore, these parameters are essential and should be presented in the manuscript. Although velocity vectors are presented in Figure 8a, it would be great to also include these vectors for all models shown in Figures 6, 7, and 9.

I would also like to invite the authors to provide videos of a representative simulation of both a permeable and impermeable sill as supplemental material. These videos could visualise fluid pressure, temperature, CH4 concentration, and velocity vectors, which would allow the reader to see how the fluid flow evolves over time and how the described phases of fluid flow form.

5) Discussion

Impact of permeable sills on hydrothermal flow

This section is a good summary of the distinct fluid flow phases observed in the models, but it misses some explanation on which mechanisms or physical parameters control each phase. For example, it is not clear, which mechanism causes a contact-parallel flow in Phase 1. In Section 4.4, the authors suggest that impermeable sills favour fluid pressure build-up and contact-parallel flow toward the sill tips. This explanation should be extended and data to support this interpretation should be included in the manuscript.

In Phase 2, pore fluid pressure (overpressure and rapid pressure dissipation after the formation of cooling joints) is used to explain the upward-directed fluid flow and hydrothermal “flushing”. As mentioned above, I would like to request the authors to provide data to support this interpretation.

It is also not clear to me how the change in pore pressure distribution initiates the vortex flow in Phase 3. Again, no data is provided to support this interpretation and the mechanism that initiates the vortex flow should be described and discussed in more detail.

As mentioned in comments 1 and 2 , permeability is a key controlling parameter that affects hydrothermal systems. This section should therefore discuss how different permeabilities would affect the described flow phases. Do the individual flow phases also form for low-permeability host rocks (k <1E-16 m²) and sill permeabilities as described by Spacapan et al. (2020)?

Implications for igneous petroleum systems in the Río Grande Valley

Based on the permeabilities used by the authors, I am not convinced that it is feasible to discuss implications for the described field location (please see comments 1 and 2).

Please find additional comments and technical corrections in the attached pdf.

Impact of permeability evolution in igneous sills on hydrothermal flow and hydrocarbon transport in volcanic sedimentary basins

The authors investigate the intrusion of magmatic sills into carbon-rich sedimentary basins that drive the production and transport of hydrocarbon-bearing fluids. They review field observations from the Neuquén Basin, Argentina, that reveal highly fractured sills containing bitumen- and calcite-bearing veins. Raman spectroscopy on the vein-filling material reveals graphitization, from which temperatures are estimated to be ~350-500 C. Based on these observations, they suggest that thermal fracturing of the sill during cooling allows the transport of pressurized hydrocarbons through the sill, which they point out is counter to the existing models that treat sills as impermeable boundaries. To investigate how the development of fractures during cooling influences the production and transport of hydrocarbons, they construct numerical models for coupled heat and fluid transport for the case of an instantaneously emplaced sill within a carbon-bearing shale host rock. 

The primary results of the modeling work are the identification of three stages of flow: 1) an early-phase contact-parallel flow regime while the sill is still hot and impermeable; 2) the development of permeability within the sill leading to vertical flow of fluids through the sill “flushing”; and 3) late-stage, slow backflow of fluids into the sill from the tip region inward. The authors show a few examples of how flow patterns differ with different sill thickness and intrusion depths constrained by the field data.

The primary conclusion is that the development of permeability in cooling sills drastically alters the flow paths and accumulation of hydrocarbons produced by heat from the sills. In a general sense this result is significant and the work should be shared with the community. However, the analysis falls short of addressing many potentially rich complexities in any kind of quantitative way. After reading this manuscript, I certainly had a greater appreciation for the number of coupled physical processes involved in these systems: cooling and fracturing of the sill, heating and maturation of the source rock, and porous and fracture-dominated fluid flow. However, I am not sure I learned a substantial amount about how these processes interact to give rise to the observed phenomena. Part of the reason for this could be the open-ended and somewhat exploratory nature of the modeling – the question of “how does permeability development influence hydrothermal flow” is too open-ended. If there were a clearly-defined question or testable hypothesis, then the modeling inputs/outputs could be analyzed in a more systematic way to address the question. For example, the authors could potentially use the model to address the following questions: Under what conditions (combinations of sill thickness, depth, initial temperatures, source rock chemistries) do we expect the development of through-going bitumen deposits in the sill? How does the accumulation of bitumen in the sill change over time and space as a function of the key parameters? (Can this then be used to compare with the field distributions?) 

Putting aside the open-endedness of the study and generally qualitative analysis of the modeling results, I still had a few questions about the treatment of permeability in the sill. I suggest that the authors address these questions prior to final publication:

Given poor constraints on the permeability of fractured media, how do assumptions about temperature-permeability relationships influence the results?

Flow through fractures is a different transport mechanism than porous flow, especially given the potential for development of permeability anisotropy related to fracture orientations. How does the treatment of fracture-dominated flow as porous flow potentially influence the flow pattern results? Would backflow be possible if fractures are vertically oriented?