| 09 Feb 2026
Placing constraints on submarine permafrost extent along the U.S. Beaufort Shelf using thermodynamic modeling
Abstract. This study constrains submarine permafrost extent along the U.S. Beaufort Shelf through eighteen thermodynamic simulations conducted along an offshore transect from Oliktok Point, Alaska, to improve characterization of the nearshore Arctic Ocean. The results demonstrate that lithology, pore water salinity, and geothermal heat flux are key controls on submarine permafrost distribution and offshore continuity. Comparison with existing geophysical observations suggests pore water salinity is near seawater levels, while geothermal heat flux could be elevated above the regional average. Synthetic resistivity logs reveal that ignoring subsurface lithologic heterogeneity can lead to misclassification of ice saturation, underscoring the need to incorporate geological context in geophysical interpretations. Moreover, comparison with prior resistivity-based interpretations suggests that some deep high-resistivity signals may reflect hydrocarbons trapped beneath the permafrost wedge rather than additional ice-bearing permafrost. Temperature profiles extracted at 2 m sediment depth provide additional constraints on subsurface thermal regimes and support calibration of Distributed Temperature Sensing (DTS) data in the absence of direct absolute temperature observations. The close agreement among modeled profiles indicates that localized DTS anomalies may represent true thermal perturbations, potentially linked to methane seepage. Collectively, these findings place meaningful limits on submarine permafrost extent along the U.S. Beaufort Shelf and provide a framework for integrating thermodynamic modeling, resistivity analysis, and DTS monitoring into future investigations of Arctic permafrost dynamics and their broader climatic implications, while also suggesting that submarine permafrost in this region may extend farther than previously recognized.
In this work, the set up of present day marine permafrost under the sea bed below a seafloor cable in the U.S. Beaufort Shelf is simulated using PFLOTRAN, by considering the mass and energy transfer within the seabed since the LGM, while taking into account of sea level variations. The performed simulations are challenging and the results are interesting, paving the way to the interpretation of the DATS data acquired offshore from Oliktok Point along the seafloor cable.
Meanwhile, the presentation of the modelling approach suffers of flaws and is sometimes not complete enough. Besides, I have a few technical concerns regarding the simulations themselves, especially the lack of convergence study. Figures should also be improved. Thus I recommend a major revision of this manuscript.
General comments:
- More details needed on the simulations set up – complete presentation of the simulation domains, of the boundary conditions, etc.
- Minor problems of structures (e.g.: simulations described before simulator, Fig 2 commented after Fig 7)
- No convergence study provided. How do look the results when using twice finer discretizations (both in space and time)? By quantifying the discrepancies between the used simulations and the finer ones, one can estimate the truncation errors, and then make sure that these truncation errors do not create artifacts important enough to mislead the analysis of the results.
Specific comments:
l 162: ‘2.1 model domain’ The simulation set up (discretization, etc) should be described after the presentation of the simulation methods, because the first ones depend on the second ones.
l 163-168: “Each geological model is discretized into ∼560,000 grid cells on an irregular, Cartesian grid. Along the horizontal X-dimension, grid cell width is 20 m. In the vertical Z-dimension, the resolution varies, ranging from 1 m near the surface of the domain to 20 m at the bottom” The convergence study performed for establishing this grid should be briefly summarized here.
l 183: Figure 4 shows three main classes of grid cell size. Why three ? How do the depths of transition between these classes have been chosen?
l 196: How the values given in Table 1 have been chosen? Here the measurements or bibliographical used for grounding the choices of these values should be given. Besides, porosity of sand varies over 10% and porosity of clay varies over 60%, but in both cases permeability is kept constant, which looks contradictory. Why not using also depth-dependent permeability?
L 207-211: Giving the initial conditions before the boundary conditions is confusing. For instance it makes difficult the understanding of the spin up. By the way, what means exactly a ‘short’ spin up?
L 216: First mention of Figure 2, after the comment of Figure 7. Figure 2 must be inserted here, or it must be commented earlier.
L 216-217: “These factors are applied at the top boundary of the model, establishing transient Dirichlet boundary conditions for pressure, temperature, and mass fraction of salt.” Figure 2 does not give any information regarding mass fraction of salt.
L 227: Then why not also considering cases with the most probable value of geothermal heat flux, 66 mW.m-2?
L 229-230: “For the lateral boundaries of the model domain, hydrostatic pressure and Dirichlet boundary conditions for the mass fraction of salt are applied.” The presentation of the simulation domains and boundary conditions is weak. For instance, what is the value for the lateral Dirichlet boundary conditions for the mass fraction of salt? A figure summarizing the geometry and the full set of boundary conditions for all primary variables (temperature, pressure and mass fraction of salt) should be provided.
L 243: I guess Dα should also be indexed with respect to the component index j. If all the components are assumed to have the same diffusion coefficient, it should be specified and the errors associated to this simplifying assumption should be briefly discussed.
L 244: wrong writing, ‘qαfor α = l, g’ would be better.
L 246: What formulation is used for the relative permeabilities kαr – Brooks-Corey, van Genuchten, anything else? This information should be given here.
L 249: Considering equation (4),
1) what Sα stand for? Only sα is introduced in the text.
2) It seems that heat conduction is considered to happen only in the solid phase (‘the rock’). Of course conduction occurs in all the phase of the porous medium, solid, ice, liquid, hydrate and gas. Is this a typo or a simplifying assumption? In the latter case the associated errors should be briefly discussed.
3) the heat capacity Cp looks also only computed by taking into account the solid phase. Same remark that for heat conduction.
L 253-256 : “The phase transition between fully liquid-saturated and partially gas-saturated is handled by a conditional that checks the dissolved gas mole fraction; if its value exceeds solubility, a gas phase forms, and the primary variables switch to gas pressure, gas saturation, and temperature.” Then I understand that in partially saturated conditions, an assumption is made regarding the flow of the liquid phase, so that the two-phase flow may still be described by a single-equation model. What is this assumption? It should be reminded here.
L 259-260 : “and a parameter defined in Grenier et al. 2018.” Please say explicitly which parameter. “and the parameter X defined in …”
Figure 11: The permafrost profiles inferred from the geophysical observations should be included in this Figure. The lithological logs should also be reminded alongside each simulations, to ease the interpretation of vertical variations.
Figure 12: it seems that there is a problem of legend – no red dashed lines, but grey ones instead. More importantly, I think that here observed resistivity should be presented as well, for allowing visual comparison between observations and simulation results.
L 400-404 : “Our research suggests that this signal … which may be misinterpreted as additional permafrost.” Here the formulation should be clarified. At the beginning the presence of hydrocarbons seems hypothetical (‘attributable’), while at the end it seems demonstrated (‘indicates’). It can’t be both at the same time.
L 406-409 : “By comparing these temperature profiles with real DTS data, assessments regarding the accuracy of the DTS measurements and determinations of whether they fall within a realistic temperature range can be made.” Indeed, comparison between observed and simulated temperature would be of great interest.
L 413 : typo ‘?’
L 416 : ‘Frederick et al, 2026 (in prep)’ – references to not already published work should be avoided.
L 419-420 : “The greatest variation is observed within the first 5 km of each model (Figure 13), where the shallow water depth makes the bottom water more sensitive to seasonal air temperature changes.” I don’t understand – these are modelling results, right? At line 220, it is specified that “When sections of the transect become submerged [...] a temperature of -1 ◦ C is applied.” Then how could seasonal variations come into play?