1 Introduction

Pockmarks are subcircular depressions on the seafloor. First reported on the Scotian Shelf (King and MacLean 1970), pockmarks have since been found worldwide in large numbers and in a wide range of settings [e.g., shallow and deep waters offshore (Hovland and Judd 1988; Pilcher and Argent 2007; Dandapath et al. 2010; Sultan et al. 2014; Roy et al. 2015; Schattner et al. 2016; Chen et al. 2015, 2017; Roelofse et al. 2020), onshore (Bogoyavlensky et al. 2020; Chuvilin et al. 2020), in lakes (Wessels et al. 2010; Reusch et al. 2015; Loher et al. 2016), and as recent and relict/buried features (Gay et al. 2007; Andresen et al. 2008; Waghorn et al. 2018; Böttner et al. 2019, 2021)]. Pockmarks range from a few meters to hundreds of meters in diameter and few centimeters to tens of meters in depth (Judd and Hovland 2007), and display variability in their shape and geometry [e.g., circular, semicircular, elliptical, and ring-like in plan-view (Sun et al. 2011; Chen et al. 2018; Zhang et al. 2020) and U-, V-, and W-shaped in cross section (Benjamin et al. 2015; Gafeira et al. 2018)]. Pockmark formation is predominantly attributed to venting of subsurface fluids to the surface accompanied by erosion and removal of sediments. Fluid sources include thermogenic gas and light hydrocarbons (Hovland et al. 1998), biogenic gas (Chen et al. 2017), groundwater (Reusch et al. 2015; Andresen et al. 2021), gas hydrate decomposition (Riboulot et al. 2011; Ruffine et al. 2013; Sultan et al. 2014), and melting permafrost (Solheim and Elverhøi 1985; Roy et al. 2012, 2015; Shakhova et al. 2017). Their formation occurs over a wide range of time scales [e.g., from hundreds of days (Cathles et al. 2010; Lohrberg et al. 2020) to hundreds of thousands of years (Judd and Hovland 2007; Wenau et al. 2017)] and is controlled by a variety of underlying geological structures and fluid release mechanisms [e.g., salt structures (Masoumi et al. 2013; Roelofse et al. 2020), faults (Pilcher and Argent 2007; Maia et al. 2016; Pape et al. 2019; Roelofse et al. 2020), ridges (Eruteya et al. 2018) and buried channels (Gay et al. 2003), sea level fluctuations (Bertoni et al. 2013; Riboulot et al. 2013), changes in sedimentation patterns, and organic matter degradation (Andresen et al. 2008; Marsset et al. 2018)]. Other processes thought to be responsible for pockmark formation include bottom current scour, reduced sedimentation, gravity flows, up-drifting ice detached from the seafloor, and biotic activity (Paull et al. 1999; Picard et al. 2018; Paull et al. 2002; Mueller 2015).

Pockmarks typically form on top of focused fluid conduits that appear as pipes or chimneys in seismic data (Andresen 2012; Karstens and Berndt 2015). Such features are generally attributed to a localized release of overpressure in the subsurface through hydraulic connection with deeper sediment layers (Hustoft et al. 2009; Cathles et al. 2010; Böttner et al. 2019). Focused fluid conduits are extremely efficient pathways for fluid migration from deeper sediments to the seafloor, and are critical for constraining global carbon emissions and predicting climate change impacts (Berndt 2005; Judd and Hovland 2007; Etiope et al. 2008; Karstens and Berndt 2015; Lohrberg et al. 2020). Moreover, fluid discharge at pockmarks leads to enhanced microbial activity, creating important marine ecological hotspots (Berndt 2005; Gay et al. 2006; Judd and Hovland 2007; Panieri et al. 2017). The presence of active gas emissions from pockmarks are reliable indicators of subsurface hydrocarbon reservoirs, and are therefore, very important for oil and gas exploration (Heggland 1998; Judd and Hovland 2007; Strozyk et al. 2018). Furthermore, actively venting pockmarks can affect seafloor stability, posing a threat to offshore infrastructure (Judd and Hovland 2007; Vanneste et al. 2014). There is also a risk of reactivation of inactive or dormant pockmarks (Judd and Hovland 2007; Pilcher and Argent 2007; Cathles et al. 2010; Roelofse et al. 2020), and consequent impacts on seafloor and subsurface installations, including subsurface carbon capture and storage (Karstens and Berndt 2015). It is, therefore, crucial to understand the processes that control the formation, progression, and activity of focused flow pathways and pockmarks, as they can provide critical proxies for subsurface fluid flow and help in the management of subsurface resources and geohazards.

From a modelling perspective, the dynamics of focused fluid flow is a highly multiphysics problem characterized by a combination of coupled and competing physical processes with many plausible driving mechanisms. In the literature, most field observations of pockmarks, pipes, and chimneys have been qualitatively explained using a variety of conceptual models based on mechanisms like capillary invasion, hydraulic fracturing, and seepage-induced erosion and weathering; or, where reactive media is involved, secondary mechanisms like local volume loss due to carbonate dissolution and hydrate dissociation are also invoked (e.g., Böttner et al. 2019, 2021; Coughlan et al. 2021; Andresen et al. 2021; Karstens and Berndt 2015; Cartwright and Santamarina 2015 and references therein).

To gain some sense of the multiphysics complexity and hierarchy of these competing mechanisms, consider a classical two-phase setting where free gas invades a fully saturated medium through a point source. In a rigid medium, the invading gas can migrate through a combination of its buoyancy and capillary forces (Helmig 1997). Capillary invasion is diffusive in nature, meaning that the capillary forces act equally in all directions around the point source, whereas the buoyancy forces direct the flow in the upward direction. If the gas source is overpressured, the nature of Darcy velocity is also diffusive, meaning that the overpressure will release equally in all directions around the point source. Therefore, in a homogeneous medium with no pre-existing preferential flow pathways (like open fractures or high-permeability lenses), buoyancy must be many orders of magnitude higher than the competing diffusive fluxes in order to focus fluid flow into chimney-like configurations (Selker et al. 2007; Cathles et al. 2010). However, if the pressure within the gas reservoir is high enough, the invading gas can initiate a sub-vertical fracture (or a network of fractures), and hydraulic fracturing can drive focused fluid flow, even for lower-buoyancy fluids (Wangen 2020; Karstens and Berndt 2015; Berndt 2005). Alternatively, the concept of solitary porosity waves has been proposed, where self-propagating, high-porosity, and high-permeability channels may emerge spontaneously due to complex and highly nonlinear coupling between fluid buoyancy, asymmetric compaction/decompaction of the pores, and viscoplastic deformations of the sediment matrix (Yarushina et al. 2015; Räss et al. 2014, 2019; Yarushina et al. 2021; Peshkov et al. 2021).

While hydro-geomechanical coupling and fracture propagation are considered the most common drivers of focused fluid flow, especially from deeper source rocks, there are limits imposed by the large overpressure magnitudes necessary for keeping the fractures open. Moreover, the upper sediment layers are typically not fully consolidated, and the mechanical energy rapidly dissipates in these layers, limiting further upward propagation. In fact, in these cases where the fracture does not propagate up to the seafloor, the tip of the fracture can itself be seen as a point gas source where the link between focused fluid flow through the upper sediment layers and emergence of pockmarks cannot be explained through hydro-fracturing or buoyancy and capillary forces alone.

Especially in such cases, the mechanism of seepage-induced erosion and hydrodynamic sediment transport (i.e., erosive fluidization) may play an important role in the genesis of focused flow pathways (Cartwright and Santamarina 2015). For example, field observations and core samples show evidence of fluidization and mobilization of sediments within focused flow conduits (Kang et al. 2016; Hurst and Cartwright 2007; Huuse et al. 2005; Böttner et al. 2021). Sandbox experiments for simulating piercement structures with air/gas injection (McCallum 1985; Nichols et al. 1994; Frey et al. 2009; Nermoen et al. 2010; Philippe and Badiane 2013) also show that seepage erosion can form surface and subsurface morphologies that are remarkably similar to pipes, chimneys, and pockmarks. Recently, Vidal and Gay (2022) have elucidated the merits and challenges of modelling the dynamics of fluid migration in sedimentary layers in the upper 500 m using concepts of frictional multiphase flow in dense deformable media (Juanes et al. 2020). In the same review, the authors have also briefly proposed conceptualization of a particles-based (i.e., discrete element method coupled with computational fluid dynamics, DEM-CFD) modelling approach.

Given the level of conceptual complexity and the hierarchy of competing multi-physics mechanisms, in this manuscript we focus our attention narrowly on the dynamics of erosive fluidization and its role in the genesis of focused flow pathways and fluid escape structures like pockmarks. More specifically, we study if, and under what conditions, this mechanism can be a primary driver of flow localization into focused flow configurations. As discussed above, in the literature there is a lack of mathematical formalism for testing this mechanism. While the model concept proposed by Vidal and Gay (2022) is promising for analyzing multi-physics couplings at small scales (e.g., sandbox), the computational costs associated with particles-based numerical schemes are formidable for realistic simulations. Therefore, we propose a mathematical model framework on the continuum scale that can resolve the flow localization through internal erosion and fluidization, and simulate the evolution of seafloor and subsurface morphology through sediment transport and redistribution. Within the scope of this study, we focus narrowly on the mechanism of erosive fluidization. However, the model is developed in a generalized manner such that it can be modularly extended to include geomechanical and chemomechanical coupling conditions. We have developed our numerical simulator in-house using an open source framework, and in this manuscript we provide a detailed description of the model as well as the simulator development. Using numerical simulations of an idealized scenario of gas escape from an overpressured shallow gas reservoir, we show how the process of erosive fluidization can localize the fluid flow into well-defined chimney-like configurations and impact the subsurface and seafloor morphologies. We also analyze the influences of the sediment properties related to flow and erosion on the morphology of the emerging fluid escape structures.

2 Methods

2.1 Model Concept

To model the physics of erosive fluidization, we conceptualize the subsurface sediment as an additive decomposition of two distinct physical states (or phases, in the macroscopic sense): (1) intact sediment, where the porous structure is preserved, and (2) fluidized sediment, where the porous structure is destroyed, and the granular material is suspended in water in a dense muddy slurry. This conceptual model is similar to the sand production models (e.g., Vardoulakis et al. 1996; Papamichos and Vardoulakis 2005; Steeb et al. 2007). Here, phase transition of the intact sediment into fluidized sediment and vice versa is controlled by erosion due to fluid seepage, and deposition due to limited carrying capacity of pore water. Therefore, fluid flow through the sediment column drives phase transitions and leads to a mass-conservative redistribution of the granular material, opening up preferential flow pathways and fluid escape structures. Sediment redistribution manifests as a variety of distinct morphological features at the seafloor. Constraining the characteristics of the internal sediment–fluid interactions can offer important insights to link seafloor morphological observations with subsurface fluid flow.

2.2 Mathematical Model

The above conceptual model is formalized through a generalized mathematical and numerical framework where the coupled fluid flow, sediment–fluid interactions, and sediment transport and redistribution are described within the macroscopic theory of porous media, and the changing seafloor morphology is resolved as a manifestation of the redistribution of aggregate sediment mass (i.e., the sum of the intact and the fluidized sediment phases).

2.2.1 Preliminaries

Our domain of interest \(\Omega \subset {\mathbb {R}}^d\) with \(d=\{2,3\}\) is partitioned into two distinct sub-domains, viz., the free surface-water domain \(\Omega _w \subset \Omega \) and porous subsurface domain \(\Omega _p \subset \Omega \), s.t. \(\Omega _w \cup \Omega _p = \Omega \) and \(\Omega _w \cap \Omega _p = \emptyset \). Furthermore, the inner boundaries between these sub-domains, denoted with \(\Gamma _{wp}\subset {\mathbb {R}}^{d-1}\), are not stationary but rather, may evolve over time due to processes such as erosion, sedimentation, surface-water runoff, and other external factors such as sea level changes. Mathematical description of conservation laws in both sub-domains is based on the homogenized variables defined over a representative elementary volume (REV). On the REV scale, the following homogenized variables are defined: local porosity \(\phi \left( \textbf{x},t\right) :=\dfrac{V_p}{V_\textrm{REV}}\) and local wetting phase saturation \(S_w\left( \textbf{x},t\right) :=\dfrac{V_w}{V_p}\), where \(V_\textrm{REV}\subset \Omega \) is the volume of an arbitrary REV, \(V_p\subset V_\textrm{REV}\) is the volume of void space where fluid flow is possible, \(V_w\subset V_{p}\) is the volume of the wetting fluid, \(\textbf{x}\in \Omega \) is the position, and \(t\subset {\mathbb {R}}\) is the time. Furthermore, the void spaces are fully occupied by at least one or both of the following two fluids: seawater and invading fluid (e.g., gas), s.t. \(S_n:=\dfrac{V_\textrm{REV}-V_w}{V_{p}}=1-S_w\). Finally, for \(\phi =0\), the definitions of \(S_w\) and \(S_n\) break down, and the model degenerates. Therefore, this model is only valid for porosity strictly greater that 0. The porous domain is characterized by \(0<\phi <1\), and the surface water domain by \(\phi =1\).

2.2.2 Conservation Equations

We consider a multiphase system composed of three distinct phases: (1) seawater, denoted with subscript w, (2) an “invading” phase (e.g., gas, light hydrocarbons) denoted with subscript n, and (3) a continuum sediment phase denoted with subscript s. The phase-wise mass balance can be summarized as

$$\begin{aligned} \partial _t \phi \rho _w S_w + \nabla \cdot \rho _w \textbf{v}_{w}&= Q_w, \end{aligned}$$
(1)
$$\begin{aligned} \partial _t \phi \rho _n S_n + \nabla \cdot \rho _n \textbf{v}_{n}&= Q_n, \end{aligned}$$
(2)
$$\begin{aligned} \partial _t \left( 1-\phi \right) \rho _s&= Q_s, \end{aligned}$$
(3)

where \(\rho _{\left( \cdot \right) }\left( \textbf{x},t\right) \) is the density, \(Q_{\left( \cdot \right) }\left( \textbf{x},t\right) \) the volumetric source terms, and \(\textbf{v}_{\left( \cdot \right) }\left( \textbf{x},t\right) \) the phase velocity for the respective homogenized phases.

Finally, we assume that the eroded sediment particles form a suspension within the pore fluids, and express the transport and eventual downstream deposition of the fluidized sediment, denoted with subscript f, as

$$\begin{aligned} \partial _t \phi \rho _w \left( S_w + \frac{\rho _n}{\rho _w} S_n\right) \Theta _{f} + \nabla \cdot \rho _w \left( \textbf{v}_{w} + \frac{\rho _n}{\rho _w} \textbf{v}_n \right) \Theta _{f}&= - Q_s, \end{aligned}$$
(4)

where \(\Theta _{f}\left( \textbf{x},t\right) \) is the mass fraction of the sediment suspended in the pore fluids. A main reason to treat sediment particles as suspended in water, and not a separate fluidized phase, is that the phase boundaries between a continuous fluidized phase and sediment-free pore water are poorly constrained.

The carrying capacity of pore fluids diminishes rapidly with decreasing phase density; for the gaseous invading phase, where \(\rho _n/\rho _w\ll 1\), we can reasonably assume that the sediments form a suspension only in the water phase, and rewrite the transport equation for fluidized sediment as

$$\begin{aligned} \partial _t \phi \rho _w S_w \overline{\Theta }_{f} + \nabla \cdot \rho _w \textbf{v}_{w} \overline{\Theta }_{f}&= - Q_s, \end{aligned}$$
(5)

where \(\overline{\Theta }_{f}\left( \textbf{x},t\right) \) is the mass fraction of the sediment suspended in seawater.

Note that Eq. 4 is a general statement of mass balance for fluidized sediment in a two-fluid-phase setting, whereas Eq. 5 is a special case applicable only to the scenarios with high density contrast between pore water and the invading phase. Within the scope of this study, we will limit analysis to gas as an invading phase, and will therefore focus on the model with Eq. 5.

2.2.3 Phase Velocities

The erosive action of subsurface fluid seepage is generally characterized by lower Reynolds numbers and larger time scales. We, therefore, assume that within this study, the Darcy model is sufficient to resolve the fluid phase velocities in both sub-domains, s.t. velocity fields evolve based only on the body forces and gradient of the respective phase pressures, \(P_\alpha \left( \textbf{x},t\right) \) for each \(\alpha =\{n,w\}\).

$$\begin{aligned} \text {For each } \alpha =\{w,n\},\quad&\textbf{v}_\alpha = - \textbf{K}\left( \phi \right) \frac{k_{r\alpha }\left( S_w\right) }{\mu _\alpha } \left( \nabla P_\alpha + \rho _\alpha \textbf{g} \right) , \end{aligned}$$
(6)
$$\begin{aligned} \text {where } \qquad&K_{i} = K_{0,i} \exp {\left[ a_{0,i} \left( \frac{\phi -\phi _0}{1-\phi _0}\right) \right] } \qquad \text {for each } i\in d,\end{aligned}$$
(7)
$$\begin{aligned} \text {and } \qquad&k_{rw}= S_w^{\left( 2/\lambda +3\right) }, \nonumber \\&k_{rn}= \left( 1-S_w\right) ^2 \left( 1-S_w^{\left( 2/\lambda +1\right) } \right) . \end{aligned}$$
(8)

The intrinsic permeability of the sediment, \(\textbf{K}\), is a d-dimensional second-order diagonal matrix. \(K_{0,i}\) and \(\phi _0\) are “reference” permeability and porosity of the intact (or uneroded) sediment, and \(a_{0,i}\) is a model parameter that controls the range of permeability variation with respect to porosity (Hommel et al. 2018). Note that when \(\phi =\phi _0\), \(K_i=K_{i,0}\), and when \(\phi =1\), \(K_i:=K_{i,\textrm{max}}=K_{0,i}\exp {\left( a_{0,i}\right) }\). Finally, \(k_{rw}\) and \(k_{rn}\) are the relative permeabilities of the w- and n-phases, respectively (Helmig 1997).

Note that, in general, a combination of Darcy–Brinkmann–Stokes models would be more accurate to resolve the flow velocities in the vicinity of the seafloor, especially when the erosion due to surface water runoff is more dominant compared to the erosive action of internal fluid seepage. However, the computational overhead of these models is extremely steep, and to the best of our knowledge, there is very limited literature on the extension of these models to multiphase settings with evolving internal sub-domain boundaries. Therefore, we do not consider a Darcy–Brinkmann–Stokes formulation within the scope of this manuscript, but acknowledge this as an important extension for more quantitative analyses in future.

2.2.4 Capillary Pressure

Capillary forces lead to a pressure jump across the fluid–fluid phase boundaries, s.t.

$$\begin{aligned} P_n-P_w := p_c\left( \phi ,S_w\right) , \end{aligned}$$
(9)

where \(p_c\) is called the capillary pressure. Here, we choose a standard Brooks–Corey (Helmig 1997) model to parameterize the capillary pressure function, and extend it using a power-law model to account for the effects of changing porosity (Goda and Behrenbruch 2011)

$$\begin{aligned} p_c = p_{e,0} \left( \frac{1-\phi }{1-\phi _0}\right) ^\beta S_w^{-1/\lambda }. \end{aligned}$$
(10)

The exponents \(\lambda \ge 1\) and \(\beta \ge 1\) are material parameters related to the particle size distribution, and \(p_{e,0}\ge 0\) is the capillary entry pressure, which is a material parameter representative of the average pore sizes in the sediment matrix.

Note that while our general model resolves capillary pressure effects through Eq. 9, within the scope of this study, \(p_c\) is ignored so that the focus of our our analysis remains only on the erosive fluidization mechanism, isolated from any consequences related to the capillary pressure hypothesis (Cathles et al. 2010). We consider this as an important step because capillary-driven flows and erosive fluidization are competing hypotheses, and since capillary pressure introduces a very strong nonlinearity, it can become nearly impossible to analyze whether erosive fluidization can drive flow localization independently or not.

2.2.5 Source Terms

In this study, there are no external sources and sinks present for the pore fluids. Furthermore, since we assume that the pore fluids are immiscible, there are also no internal sources and sinks due to phase transitions and mass exchange. Therefore, \(Q_w=0\) and \(Q_n=0\).

However, due to the erosion and deposition processes, there is a continuous exchange of mass between intact sediment and fluidized sediment, which leads to the internal source terms

$$\begin{aligned} Q_s := \sum \limits _{\alpha =n,w} \left( -\epsilon _{\alpha }\left( \phi ,|\textbf{v}_\alpha |\right) + \delta _{\alpha }\left( \phi S_\alpha ,\Theta _f\right) \right) , \end{aligned}$$
(11)

where \(\epsilon _\alpha \) and \(\delta _{\alpha }\) are phase-wise erosion and deposition rates, respectively (see Rahmati et al. 2013 and references therein)

$$\begin{aligned}&\epsilon _\alpha = e_{\alpha ,0} \left( 1-\phi \right) ^m |\textbf{v}_\alpha |^n, \end{aligned}$$
(12)
$$\begin{aligned}&\delta _\alpha = ds_{\alpha ,0} \left( \frac{\Theta _f}{\phi S_\alpha }\right) ^\gamma , \end{aligned}$$
(13)

internal erosion length scales \(e_{\alpha ,0}\), intrinsic deposition rates \(ds_{\alpha ,0}\), and empirical parameters m, n, and \(\gamma \).

As previously stated, in this manuscript we focus only on the invasion of the gas phase, s.t. \(\frac{\rho _n}{\rho _w} \ll 1\). Therefore, the general source terms in Eqs. 12 and 13 reduce to the special cases

$$\begin{aligned}&\epsilon _n = e_{n,0} \left( 1-\phi \right) ^m, \quad \epsilon _w = 0, \quad \delta _w = ds_{w,0} \left( \frac{\overline{\Theta }_f}{\phi S_w}\right) ^\gamma \quad \text {and} \quad \delta _n = 0. \end{aligned}$$
(14)

2.3 Numerical Solution Scheme

The mathematical model has four main governing equations: (1)–(3) and (4 or 5). We chose the following primary variables for our numerical model: \(P_w\), \(S_n\), \(\phi \), and \(\Theta _f\text {(or }\overline{\Theta }_f)\).

The governing equations are discretized spatially using a fully upwinded cell-centered finite volume scheme with a linear two-point flux approximation defined on an orthogonal mesh with rectangular \((d=2)\) or cuboidal \((d=3)\) elements, and temporally using an implicit finite difference scheme. Additionally, in order to evaluate the erosion rates, the fluid velocity fields are reconstructed using an L2 projection (http://www.mathematik.tu-dortmund.de/~featflow/en/software/featflow2/tutorial/tutorial_l2proj.html) of the phase pressures from their native P0 space to a higher Q1 space. The spatially discretized model is partitioned into three sub-modules: (1) the two-phase flow module \({\mathcal {M}}_1\), composed of the governing equations (1) and (2) with primary variables \({\mathcal {P}}_1=\left[ P_w,S_n\right] ^T\), (2) the L2 projection module \({\mathcal {M}}_2\), with the projected phase pressures as “intermediate” primary variables, \({\mathcal {P}}_2=\left[ \overline{P}_w,\overline{P}_n\right] ^T\), and (3) the sediment module \({\mathcal {M}}_3\), composed of the governing equations (3) and (4 or 5) with primary variables \({\mathcal {P}}_3=\left[ \phi ,\Theta _f(\text {or }\overline{\Theta }_f)\right] ^T\). At each time step, the numerical solution of the coupled problem is obtained by solving the sub-modules iteratively in a blocked Gauss–Seidel scheme (Gupta et al. 2015).

The resulting numerical scheme is implemented within version 2.8 of the DUNE-PDELab framework based on C++ (Bastian et al. 2010; Sander 2020) and uses the in-built matrix assembler, linearization algorithm (Newton method with a numerical Jacobian), and linear solver [parallel algebraic multi-grid (AMG) solver with a stabilized bi-conjugate gradient (bi-CG) preconditioner]. The DUNE libraries used in this study are preserved at https://gitlab.dune-project.org/pdelab/dune-pdelab and developed openly at https://www.dune-project.org/. The source code for the model and test scenarios presented in this manuscript is archived in the following public repository: https://github.com/shub-G/ErosiveFluidizationModel/releases/tag/v1.0.0. The computations for this study were performed on the high-performance computing cluster at Kiel University (CAU).

2.4 Test Setting and Computational Domain

Fig. 1
figure 1

a Test setting and two-dimensional computational domain, b the representative control volume and definitions of the homogenized variables, and c model parameters used for the numerical simulations

Full size image

To understand the effects of erosive fluidization on the sediment and seafloor morphology, we simulated and analyzed the sediment redistribution in a representative geological scenario that is commonly linked with focused fluid flow and escape structures, namely, the infiltration of buried overpressured gas into an overlying sediment layer connected to the seafloor (Judd and Hovland 2007; Cathles et al. 2010; Cartwright and Santamarina 2015). The idealized test setting, shown in Fig. 1, considers a light hydrocarbon (e.g., methane gas) trapped under high pressure in a buried reservoir sealed by a capillary barrier. The overlying sediment is assumed to be stratigraphically homogeneous, fully water-saturated, and continuously connected to the seafloor. At \(t=0\), a fracture spontaneously punctures the capillary barrier and allows the overpressured gas to escape. Note that this fracture does not cross into the porous sediment layer. Rather, it only opens a flow pathway for overpressured gas into the overlying sediment layers. The overpressure in the gas reservoir is not maintained indefinitely, but only until the gas plume reaches the seafloor. We identify a two-dimensional computational domain \(\Omega \) as a region around the fracture opening above the capillary barrier. The computational domain encompasses the overburden \(\Omega _{p}\) and the surface water column \(\Omega _w\), and explicitly resolves the seafloor \(\Gamma _{wp}\). The computational domain parameters and the material properties are summarized in Fig. 1. The reference permeability and porosity of the sediment are chosen to represent permeable silty sands (Forster et al. 2003).

It is important to note that we do not resolve the source of the free gas in the gas reservoir and the cause of the fracture. We also ignore the effects of bottom water currents in \(\Omega _w\) in the vicinity of \(\Gamma _{wp}\). These assumptions allow us to focus exclusively on the erosive fluidization and flow localization processes within \(\Omega _p\), and isolate the evolution of \(\Gamma _{wp}\) from drivers other than the erosion in the subsurface. Similarly, we also ignore the capillary effects and geomechanical feedback so that we can categorically identify the correlation between the parameters of erosive fluidization and the emerging subsurface and seafloor morphologies.

3 Results and Discussion

Fig. 2
figure 2

Impacts of erodibility \(r_0 \in \left\{ 100, 50, 25\right\} \) on a total sediment distribution \(\left( s:= \left( 1-\phi \right) + \frac{\rho _w}{\rho _s}\phi S_w \overline{\Theta }_f\right) \), b sediment distribution in pore water \(\left( s^\prime := \frac{s}{S_w} \right) \), c gas saturation \(\left( S_n:=1-S_w\right) \), and d gas pulses on the seafloor. The numerical results are plotted for the particular scenario with flow anisotropy \(K_F=1\) and correspond to time \(t=2\) years. The quantity s shows the gas chimney/pipe encased in a low-porosity, low-permeability halo in the subsurface, culminating in a W-shaped pockmark on the seafloor. The quantity \(s\prime \) highlights the gas flow pathways, which appear to be annular within the pipe/chimney. High \(r_0\) leads to a narrower focused flow path with a tighter sediment halo (i.e., pipe), while low \(r_0\) leads to a wider focused flow path with a more diffuse sediment halo (i.e., chimney). Moreover, a combination of low anisotropy (\(K_F=1\)) and high erodibility (\(r_0=100\)) leads to a ring-shaped pockmark. In general, active fluid flow tends to form W-shaped pockmarks. Finally, sediment erodibility also impacts gas flow, where high \(r_0\) leads to high-frequency and low-amplitude gas pulses

Full size image

Selected results from the numerical study are shown in Figs. 2 and 3 where snapshots of (a) redistributed sediment (i.e., volume of fluidized and intact sediment per unit REV, \(s:=\left( 1-\phi \right) +\frac{\rho _s}{\rho _w} \phi S_w \overline{\Theta }_f \)), (b) redistributed sediment relative to pore water (i.e., volume of fluidized and intact sediment per unit pore water, \(s\prime :=\frac{s}{S_w} \)), and (c) gas saturation \(S_n\) are plotted at \(t=2\) years, and (d) average \(S_n\) measured at the seafloor (black) and \(S_n\) at the deepest point of the pockmarks on the seafloor (green) are plotted over the duration of the simulations. The quantity of interest (qoi) s is a conservative property bounded between 0 and 1. It shows the distribution of the aggregate sediment mass (i.e., sum of intact and fluidized sediment) and highlights the focused flow pathways. The qoi \(s\prime \) is not a bounded quantity, but it is very interesting because it highlights the gas channels, which cannot be identified by looking at the distribution of aggregate sediment mass alone. Results show that the mechanism of erosive fluidization leads to localization of gas flow into distinct focused flow pathways in the subsurface, which manifest as pockmarks at the seafloor, as marked in Figs. 2a and 3a. The shape and size of the focused flow pathways and pockmarks are analyzed with respect to three main sediment characteristics related to seepage and erosion: (1) erodability \(r_0:=\frac{e_{n,0}}{ds_{w,0}}\) (i.e., ratio of erosion and deposition rates; see Fig. 2), (2) flow anisotropy \(K_F:=\frac{K_{0,0}}{K_{0,1}}\) (i.e., lateral vs. vertical permeability; see Fig. 3), and (3) intrinsic permeability \(K_{0,1}\) (see Fig. 5).

Fig. 3
figure 3

Impacts of flow anisotropy \(K_F \in \left\{ 1, 10, 100\right\} \) on a distribution of total sediment mass \(\left( s:= \left( 1-\phi \right) + \frac{\rho _w}{\rho _s}\phi S_w \overline{\Theta }_f\right) \), b distribution sediment mass in unit pore water volume \(\left( s^\prime := \frac{s}{S_w} \right) \), c gas saturation \(\left( S_n:=1-S_w\right) \), and d gas pulses on the seafloor. The numerical results are plotted for the particular scenario with sediment erodibility \(r_0=50\), and correspond to time \(t=2\) years. The quantity s shows the gas chimney/pipe encased in a low-porosity, low-permeability halo in the subsurface, culminating in a W-shaped pockmark on the seafloor. The quantity \(s^\prime \) highlights the gas flow pathways, which appear to be annular within the pipe/chimney. High flow anisotropy makes the sediment halo more diffuse due to higher lateral sediment transport. The impact of flow anisotropy on pockmark shape is nonlinear because of a competition between lateral sediment transport and limited carrying capacity of the pore water. Flow anisotropy has a large impact on the gas flow, where high \(K_F\) leads to high-frequency and high-amplitude gas pulses

Full size image

The key findings based on our numerical results are as follows:

  1. 1.

    Erosive fluidization leads to morphological features like a conical focused-flow pathway with annular gas flow, which may be interpreted as gas pipe or chimney, encased in a halo of low permeability sediment that acts as an effective seal against lateral gas transport. The sediment within the focused-flow pathway undergoes intense fluid seepage-driven mixing, which results in brecciation (interpreted in our model as regions with low intact sediment volume and high fluidized sediment volume). Figures 2a and 3a show the focused-flow pathway and the sediment halo, while Figs. 2b and 3b highlight the annular gas channels embedded within the brecciation zones. Interestingly, each of these features bears striking similarities with the sandbox experiments where the formation of piercement structures was analyzed through controlled fluidization using pressurized air injection (McCallum 1985; Nermoen et al. 2010).

  2. 2.

    Erosion and deposition are competing processes with a complex feedback loop: within the sediment, high erodability \(\left( r_0\right) \) leads to a more prominent near-cylindrical gas chimney and “tight” halo around the focused-flow pathway. On the other hand, high flow anisotropy \(\left( K_F\right) \) leads to a wider focused-flow pathway and wider and a more diffuse sediment halo.

    In the literature, the terms “chimney” and “pipe” are often used interchangeably to describe focused-flow pathways, although some authors (e.g., Andresen et al. 2021; Karstens and Berndt 2015) consider a stricter nomenclature where pipe refers to cylindrical flow conduits with sharp boundaries between the focused flow zone and the the host sediment, and chimney refers to irregular conduits with chaotic transition towards the host sediment. Based on this nomenclature, our results indicate that high \(r_0\) leads to pipes, and high \(K_F\) leads to chimneys, although the transition between the two structures is a continuous spectrum.

    On the seafloor, high \(r_0\) leads to narrower pockmarks with sharp depressions, while \(K_F\) shows more nonlinear trends related to depth and width of pockmarks. In general, high \(K_F\) implies a higher lateral flow of the sediment mass, and vice versa. If erosion rate is high, sediment tends to collect on the rim of the pockmark, forming a raised ring-shaped encasing. On the other hand, when erosion rates are low, less sediment mass reaches the surface, and therefore, ringed pockmarks do not form. For high erosion rates, if anisotropy is low, the ring shape of the pockmark is more pronounced. As anisotropy increases, more lateral flow occurs, leading to flattening of the ring, formation of secondary rings, and widening of the pockmark depressions. For low erosion rates, increasing anisotropy leads to widening of the pockmark depression up to the point where the carrying capacity of water is unable to keep up with the erosion process, after which a further increase in sediment anisotropy leads to progressively narrower and shallower pockmark depressions.

  3. 3.

    Gas release occurs in pulses (see Figs. 2d and 3d). Low \(r_0\) and high \(K_F\) values lead to high amplitude of the gas pulses (black curves), while high \(r_0\) and high \(K_F\) values result in high frequency of gas pulses. The cause of this pulsed gas flow is the nonlinear coupling between seepage velocity, erosion and deposition rates, and porosity: The rate of erosion is proportional to seepage velocity, which is proportional to porosity. Higher seepage velocity leads to higher erosion, which leads to increased porosity and even higher seepage velocity. However, higher erosion also leads to higher concentration of suspended granular material, which increases the fluid density as well as deposition rates, thereby reducing the porosity and, in turn, the seepage velocity. Since the mass transport and the sediment phase transitions are dynamic and rate-based (i.e., not spontaneous), the seepage velocity shows a periodic increase and decrease, leading to pulsed gas flow. In nature, the periodicity of higher-frequency gas pulses is likely masked by the bottom water currents. However, our results suggest that the pulsed release of gas is an intrinsic feature of the physics of flow localization through erosive fluidization. Interestingly, this pulsed gas release was also reported in sandbox experiments of piercement structures through air injection (Nermoen et al. 2010), and the localization of gas in the subsurface predicted by our model (e.g., Fig. 3c with \(K_F=100\)) closely resembles the air accent imaged in these experiments. Pulsed gas release was also observed at the Scanner pockmark (Callow et al. 2021). While at the Scanner site, gas pulses have been linked to tidal fluctuations (Li et al. 2020), our results suggest that gas pulses can, in theory, also occur due to intrinsic dynamics of flow localization without any external forcings like storm or tidal waves.

  4. 4.

    Differences in sediment–fluid interaction characteristics lead to diversity in pockmark shapes and sizes. Within the constraints of this test setting (i.e., short time scale and continuous gas seepage), two interesting pockmark geometries emerge: W-shaped and ring-shaped, as shown in Fig. 4. Although our test setting was theoretical, there are many real-world analogs where our simulation results seem to apply.

    1. 4.1.

      W-shaped pockmarks (Lazar et al. 2019; Hovland and Judd 1988; Callow et al. 2021; Watson et al. 2020; Gafeira et al. 2018), also referred as “inverted” pockmarks (Lazar et al. 2019), are associated with active fluid escape (Lazar et al. 2019; Schattner et al. 2012).

      Our results suggest that the depression of the annular gas channel follows radial symmetry.

      A real-world analog is a radially symmetric W-shaped pockmark reported in the continental shelf offshore northern Israel (Schattner et al. 2012), formed in a geological setting strikingly similar to the idealized setting considered in this study. The pockmark is roughly 60 m across and 10 m deep, and is the result of active venting of methane gas, with seabed gas emission occurring at shallow water depths (\(<100\) m). It lies in a region of low cohesion and is located directly above a chimney that extends up to the last glacial maximum (LGM) unconformity at depth between 100 and 200 m (Schattner et al. 2012).

    2. 4.2.

      Ring-shaped pockmarks have one or more concentric raised rims, formed due to a combination of high erosion rate and low flow anisotropy, as discussed in key finding no. 2. A prominent example of ring-like pockmarks is found in the Hudson Bay (Roger et al. 2011 and references therein). Movement of icebergs may have caused the breaching of the capillary seals, allowing the escape of hitherto unknown hydrocarbon fluids from the source rocks supposedly lying between depths of 80 m upto 200 m (Zhang 2008 and references therein).

      The origin of these ring pockmarks remains speculative, although similar ring-like features have been previously observed in the Hudson Basin (Dimian et al. 1983) and attributed to hydrocarbon escape from possible salt-related, block-faulting, and chemosynthetic formation mechanisms. Regardless of the underlying fluid source, the geological scenario bears similarities with our computational setting, and our simulations suggest that erosive fluidization from the migration of the escaped fluid in sediments with low flow anisotropy and high erosion rate can lead to such ring-like pockmarks on the seafloor.

  5. 5.

    The intrinsic permeability \(K_{0,1}\) does not impact the shape and size of pockmarks and pipes/chimneys, as shown in Fig. 5a (compare left and right). It only affects the timescale of their evolution (e.g., pockmarks emerge within \(\sim 1\) yr in sediments with \(K_{0,1}=10^{-13}\,\text {m}^2\) and \(\sim 100\) years when \(K_{0,1}=10^{-15}\,\text {m}^2\)). However, differences do appear in the frequency of gas pulses, with lower \(K_{0,1}\) leading to higher frequencies. This behavior is illustrated in Fig. 5b. In the same figure (5a, compare left and center), we also see that the morphology of pockmarks and pipes/chimneys is controlled by the ratio of \(e_{n,0}\) and \(ds_{w,0}\) (i.e., \(r_0\)) but not by their individual magnitudes. However, similar to \(K_{0,1}\), the magnitudes of \(e_{n,0}\) and \(ds_{w,0}\) do affect the amplitude of the gas pulses, with lower \(e_{n,0}\) (or conversely higher \(ds_{w,0}\)) leading to higher amplitude.

Fig. 4
figure 4

Comparison of selected numerically simulated pockmarks on the basis of shapes. a W-shaped pockmarks, b ring-shaped pockmarks

Full size image
Fig. 5
figure 5

Impacts of intrinsic permeability \(K_{0,1}\) and erosion rate constants \(e_{0}\) on the evolution of surface and subsurface morphological features as well as the gas flow behavior. Snapshots for scenarios with \(K_{0,1}=10^{-15}\,\text {m}^2\) correspond to time \(t=100\) years, while those with \(K_{0,1}=10^{-13}\,\text {m}^2\) correspond to \(t=1\) year. Results of a sediment distribution show identical profiles, implying that \(K_{0,1}\) has no effect on morphology, and neither does the magnitude of \(e_{0}\) (for a given \(r_0\)). However, snapshots of b gas saturation show that while the subsurface flow of gas remains identical, \(K_{0,1}\) and \(e_{0}\) have noticeable impact on gas pulses at the seafloor

Full size image

4 Conclusions

We presented a continuum-based model for simulating the formation of focused flow and escape structures based on the mechanism of erosive fluidization. Numerical simulations of an idealized scenario of gas escape from an overpressured gas reservoir showed that erosive fluidization and sediment transport lead to (1) formation of conical focused-flow conduits with brecciated core and annular gas channels encased within a halo of low-permeability sediment, (2) pockmarks of diverse shapes and sizes on the seafloor, including W and ring shapes, and (3) pulsed release of gas. Analysis of sediment characteristics revealed that intrinsic permeability has no impact on the subsurface and seafloor morphologies. Sediment erodibility and flow anisotropy emerged as dominant controls. Although the test setting was theoretical, our results have striking real-world analogs in nature as well as controlled experiments.

An important takeaway of this study is that the mechanism of erosive fluidization can lead to flow localization into focused-flow configurations without any capillary and/or geomechanical feedback, or in other words, erosive fluidization can be a primary driver of flow localization. Therefore, quantitative models for analyzing focused fluid flow must include this process along with other chemo-hydro-geomechanical multiphysics couplings. To that effect, this manuscript proposes a flexible mathematical formalism that can be integrated into other existing models and software with relative ease.

While the focus of this study was narrowly on the mechanism of erosive fluidization, it is clear that no single mechanism can universally describe all the observed features of focused fluid flow. The quantitative analysis of this problem demands a comprehensive multiphysics model framework that can handle both hydrogeomechanical processes and sediment fluidization, and possibly also biogeochemical feedback. Therefore, the similarities of our results with real-world analogs must be viewed as solely qualitative, with the caveat that both the geomechanical and capillary effects were ignored.

The computational framework for the proposed model was developed keeping the strongly multiphysics character of the focused flow problem in mind, and is therefore highly general and modular by design. In this manuscript we have demonstrated the capability of this framework in terms of handling the sediment phase transitions and transport. Future work will include more complex couplings like capillary effects, mechanical deformation, and gas miscibility.