# Evolution model for the Absheron Mud Volcano: from stratified sediments to fluid mud generation

Arthur Blouin, Nabil Sultan, Alexandra Pierron, Patrice Imbert, Jean-paul Callot

► To cite this version:

Arthur Blouin, Nabil Sultan, Alexandra Pierron, Patrice Imbert, Jean-paul Callot. Evolution model for the Absheron Mud Volcano: from stratified sediments to fluid mud generation. Journal of Geophysical Research: Earth Surface, 2020, 125 (12), pp.e2020JF005623. 10.1029/2020JF005623. hal-04202657

HAL Id: hal-04202657

https://hal.science/hal-04202657v1

Submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Evolution model for the Absheron Mud Volcano: from stratified sediments to fluid
mud generation

Arthur Blouin¹,²,³, Nabil Sultan¹, Alexandra Pierron¹, Patrice Imbert³, Jean-Paul Callot³

¹ Ifremer, Lab. Aléas Géologiques et Dynamiques Sédimentaires, Géosciences Marines, REM, Centre de Bretagne, Pointe du Diable, 29280 Plouzané.

² R&D/EP Total S.A., Avenue Larribau, 64000 Pau.

³ Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, TOTAL, Avenue de l'Université, 64013 Pau cedex, France

Corresponding author: Arthur Blouin (a.blouin@fugro.com)

Key Points:

* Simulation of mud generation is possible through 2D basin modeling including the effect of gas

* Mud ascent through gas-expansion-driven density inversion is possible

* A semi-quantitative formation model for the Absheron mud volcano is derived from modeling results.
16 **Abstract**

17 Submarine mud volcanoes are one of the most spectacular methane expulsion features at the
18 seafloor and they represent a significant geohazard worldwide. In this work we focus on the
19 physical processes controlling the initiation and early evolution of the Absheron mud volcano
20 (AMV). Our analyses were carried out based on basin modeling calibrated thanks to existing
21 seismic interpretation of the AMV, analysis of sediment samples from seabed, and data from two
22 exploration wells. Acquired laboratory geotechnical data allowed us to derive laws considering
23 the impact of gas exsolution on host sediment behaviors.

24 In this study, we identified key geological and physical conditions that led to mud volcano
25 formation: by coupling diffusion processes with hydrofracturing and fluid advection, we were
26 able to simulate the conditions required to generate mud 3.5 km below the AMV. Mud
27 remobilization up to the seabed was reproduced by using Navier-Stokes equations modified to
28 account for the impact of gas expansion on mud density. Considering density inversion only,
29 simulations indicate that mud would be extruded at the seabed 100 years after its generation, an
30 ascent rate similar to extrusion rates measured at the active Kotyrdag mud volcano in Azerbaijan.

31

32 **Plain language Summary**

33 Mud volcanoes build-up at the Earth surface from liquidized sediments remobilized from depths
34 reaching several kilometers. They are distributed globally, onshore and offshore, in different
35 geological backgrounds and represent serious geohazards for people and infrastructures.
36 Moreover, they are also among the most spectacular methane expulsion features at the surface,
37 participating to the global budget of green-house gases emissions. Their understanding is thus
38 paramount to prevent natural catastrophes such as Lusi eruption in Indonesia in 2006 and to
39 quantify the human impact on global warming. If mud volcanoes are known for centuries (Pline
40 the Elder) and studied for decades, little is known about the processes controlling mud
41 generation from stiff deep sediments and its vertical migration towards the surface. The methane
42 is often considered as being a driving element in this geological process and bubble formation in
43 compacted sediments during laboratory testing led to strong damage of sediment properties.

44 Here we develop two independent numerical models, based on simple physical processes and
45 relying on in-situ observations of the Absheron Mud Volcano (South Caspian Basin) and
46 laboratory testing. They allow to simulate mud generation at depth and its remobilization to the
47 surface. Our research work provides reliable information to describe the initiation and early
48 evolution of an active mud volcano.
1. Introduction

While extensive work has focused for decades on mud volcano (MV) architecture and on gassy sediment structure and behavior, few studies have numerically explored MV formation and evolution. Several authors have attempted to reproduce surface and subsurface MV architecture as well as their suspected morphological evolution through analogue models (Dupuis, 2017; Mourgues et al., 2012; Nermoen et al., 2010; Odonne et al., 2020; Woolsey et al., 1975). The results provide understanding of the evolution of surface structures and morphology as well as exploring possible subsurface mechanisms responsible for the observed data. However, most of these models used non-cohesive materials to simulate fine-grained cohesive sediments (Dupuis, 2017; Nermoen et al., 2010; Woolsey et al., 1975), following the theory on scaled models (Hubbert, 1937). The processes of material remobilization in these models may be different from those applying to natural MVs. Other authors developed numerical models to explain the subsurface mechanism and simulate the remobilization of mud (Brown, 1990; Collignon et al., 2018b; Gisler, 2009; Zoporowski & Miller, 2009). Rudolph et al. (2011) proposed a model where mud remobilization allows the mud chamber to extend laterally from an initial volume of remobilized mud. These models are based on physical, chemical and mechanical laws applied to natural conditions, therefore approaching processes more accurately. However, Collignon et al. (2018a) stress the difficulties of modeling clastic eruptions, as they involve different temporal and physical scales into a multi-phase and multi-process system. Additionally, none of the cited studies tackles directly the modeling of initial mud remobilization process. Deville et al. (2010) computed regional fluid flow models in the Barbados prism and showed that MVs are located above overpressured zones. Therefore, modeling of mud volcanoes focused on reproducing the observed structures or the mud flow dynamics, while other studies looked into the regional background to explain MV location. To our knowledge, none of the previous studies looked into both the regional and local processes. The impact of gas exsolution and expansion on subsurface sediment mobilization has been inferred from geophysical observations (Brown, 1990; Imbert et al., 2014; Lafuerza et al., 2012; Lawrence & Cartwright, 2010; Riboulot et al., 2013), but the physics of gas-related mud generation commonly ignored or simplified.

The Absheron mud volcano (AMV) is located in the South Caspian Basin (SCB; Figure 1). The SCB is a complex area that has recorded intense tectonic activity since the Cretaceous and which is one of the sedimentary basins having recorded the largest subsidence and sedimentation rates worldwide, notably thanks to the inception of a subduction along the Absheron ridge and through the interaction of three paleo-river deltas (Allen et al., 2002; Brunet et al., 2007; Egan et al., 2009; Golonka, 2007; Morton et al., 2003; Smith-Rouch, 2006). It is also a province of intense O&G (oil and gas) activity. This is favored by the presence of a mature Oligocene-Miocene source rock (Maykop Formation), numerous regional reservoirs deposited during the Pliocene (e.g. NKG, Balakhany, Fasila, Sabunchi and Surakhany Formations), structural traps in anticline crests and a good seal interval composed by the entire Quaternary section (Alizadeh et al., 2017; Guliyev et al., 2011; Javanshir et al., 2015). The SCB also holds the densest distribution of MVs worldwide (Dimitrov, 2002; Kopf, 2002; Milkov, 2000). They are related to gas accumulations along anticline crests and their formation seems intimately linked to fold growth (Dupuis, 2017; Fowler et al., 2000; Yusifov & Rabinowitz, 2004).

The Absheron anticline formed due to the subduction initiated at the end of Miocene (Allen et al., 2002; Golonka, 2007; Hollingsworth et al., 2008). It is cored by a deeply-rooted thrust and by related normal faults on its extrados (Blouin et al., 2019a). As it hosts a giant gas condensate
94 field, the AMV area has been densely surveyed for hydrocarbon exploration purposes and for 95 geohazard assessment (Contet & Unterseh, 2015; Dupuis, 2017; Gautherot et al., 2015; Unterseh 96 & Contet, 2015). A dense and varied dataset was accumulated over the exploration phase, with 97 notably high resolution industrial geophysical data (Figure 1), geotechnical in situ measurements, 98 sediment cores and two exploration wells together with industrial logs, giving a dense and high- 99 quality multidisciplinary dataset. Thus, the AMV can be apprehended through different 100 approaches. Based on geophysical methods, its morphology has been interpreted in detail and 101 compared with in situ data and sediment analysis in order to get a comprehensive model of its 102 formation (Blouin et al., 2019a). The sediments have been tested in laboratory through a technique 103 developed to understand the mechanical behavior of compacted sediments hosting gas exsolution 104 and to study mud generation (Blouin et al., 2019b). Finally, the comprehensive formation model 105 and the laboratory results can be confronted through numerical models computing several physical 106 and geological processes. All in all, the AMV is located in an area that gives a unique opportunity 107 to study an active MV and model its formation and evolution within its geological background

108 Previous modeling work carried out by Blouin et al. (2019a) allowed to explain the location 109 of the AMV at the crest of the Absheron anticline, coupling sedimentation-related overpressure 110 generation and 2D-diffusion equations for fluid flow and methane diffusion over 5 million years 111 (My). It showed that the AMV was formed at the location where critical overpressure (near- 112 fracture conditions) affected methane-saturated areas. However, further structural and strata pre- 113 conditioning are needed to reach hydrofracturing conditions at the AMV location and model its 114 formation.

115 The present work assesses the following and complementary main questions:

116 • Is it possible to simulate mud generation conditions at depth by considering 117 mechanical properties of sediments, sedimentation rates, structural elements and 118 the impact of gas exsolution?

119 • Is the sole impact of gas expansion on the mud properties able to drive mud up to 120 the seafloor?

121 **2. Material and methods**

122 In addition to available data from literature, we constrained the study using information 123 from two exploration wells (pressure, stratigraphy) and subsurface geometry from the 124 interpretation of a 3D seismic volume (Blouin et al., 2019a) and using laboratory testing results 125 showing the impact of gas exsolution on the mechanical properties of compacted sediments 126 (Blouin, et al., 2019b). Two different numerical models were used to model mud generation at 127 depth and to simulate local processes related to mud remobilization towards the surface.

128 **2.1. Mud generation**

129 **2.1.1. Sediment properties and geological structural geometry**

130 The experimental method used to characterize sediment properties (including 131 compressibility and hydraulic conductivities) of the modelled stratigraphic layers was already 132 discussed in Blouin et al. (2019a) as well as the software used to calculate the sedimentation- 133 related overpressure (SeCoV3 Ifremer in-house software).

134 The geometrical model, based on work by Green et al. (2009), was slightly modified
(Figure 2) to take into account the Anhydritic Surakhany Formation (ASF). The ASF is a 300 m thick interval mainly composed of interstratified evaporites and over-pressured and under-compacted clay-rich layers that displays sealing capacities on well data (Dupuis, 2017; Javanshir et al., 2015). A fault network geometry based on seismic interpretation (Blouin et al., 2019a) was also considered.

In order to optimize the calculation time, the spatial resolution was reduced by considering a vertical resolution of 192 m with a total of 50 vertical nodes, and a horizontal resolution of 244 m with a total of 225 horizontal nodes (Figure 2).

## 2.1.2. Hydrofracturing: physical processes and working hypotheses.

The 2D-diffusion method was modified from the initial version (Blouin et al., 2019a) in order to consider two additional phenomena: hydrofracturing and gas compressibility and exsolution (Boyle's and Henry's laws). In a context where the maximum principal effective stress is vertical ($\sigma'_1 = \sigma'_v$), hydrofracturing is expected to occur when and where the fluid pressure exceeds the sum of the tensile strength of the host sediment and the minimum horizontal stress ($\sigma'_3$). This simplified condition is verified on the fold extrados where normal faulting was observed and where well data show a fracture pressure smaller than $\sigma'_v$ (Blouin et al., 2019a). Nevertheless, it may not apply away from the anticline crest, the overall tectonic regime in the SCB being compressive. For the purpose of modeling, a realistic approximation is to consider that hydrofracturing occurs when pore pressure exceeds the minimum horizontal stress which corresponds to a ratio between overpressure and vertical effective stress ($\Delta u/\sigma'_v$), generally obtained from pressure logs. When the ratio $\Delta u/\sigma'_v$ reaches this critical condition, hydrofractures are generated in the concerned area. As a working hypothesis, the hydrofractured area permeability is calculated from hydraulic conductivity laws using a void ratio of 1. Methane diffusivity in the hydrofractured area is considered as being the same as in faults. Fracturing may trigger a decrease of overpressure in the fractured zone leading to a decrease in the saturation concentration and consequently to gas exsolution. Hence, the degree of gas saturation ($S_g$) after exsolution is calculated in the fractured area. The calculation follows Boyle's law (methane compressibility; Boyle, 1965) and Henry's law (dissolution/exsolution capacity of methane; Henry, 1803), both being pressure and temperature dependent. Therefore, temperature distribution was calculated over the modeled sedimentary column. We estimated a thermal gradient from temperature measurements in the two exploration wells (16°C/km) and seafloor temperature from Diaconescu et al., (2001) at 5.9°C.

## 2.1.3. Sediment damage caused by gas exsolution

Several studies noted that during marine sediment sampling, the structure and mechanical properties of sediments are impacted by gas exsolution and expansion (DeGroot et al., 2010; Esrig & Kirby, 1977; Priest et al., 2014; Sultan et al., 2010, 2012). Sultan et al. (2012) determined that damage of sediments from the Gulf of Guinea through a decrease in preconsolidation pressure (maximum effective vertical stress that a soil sample has sustained in the past) was directly linked to the degree of gas saturation. This decrease in preconsolidation pressure with increasing degree of gas saturation was already highlighted on other sediments by several studies (Hight et al., 2002; Lunne et al., 2001). Blouin et al. (2019b) observed the same behavior on sediments from the Absheron mud volcano, mud being generated for a degree of gas saturation of 38%.
The degree of gas saturation is consequently used to calculate the impact of gas exsolution on sediment damage through the equation obtained from laboratory testing on the AMV sediments (Blouin, et al., 2019b):

$$ \frac{\sigma'_{p2}}{\sigma'_{p0}} = \exp(-0.07 \cdot S_{gmax}) \qquad (1) $$

with $\sigma'_{p2}$ the preconsolidation pressure of sediments calculated after gas exsolution, $\sigma'_{p0}$ the preconsolidation pressure of sediments before gas exsolution and $S_{gmax}$ the maximum degree of gas saturation of the tested sample.

### 2.1.4. Modelling of mud generation: limit conditions

Methane and pore-pressure diffusion calculations may be initiated at different simulation stages accounting for hydrocarbon generation from mature source rocks since the beginning of Pliocene times, at the very beginning of the modeled sedimentation process. Boundary conditions, displayed in Figure 2, are defined as follows:

* Conditions 1-2: no lateral or vertical exchange of pore pressure and dissolved methane outside of the model limits.

* Condition 3: sea water has a fixed pore pressure of 0 kPa and a methane concentration of 1.0 x 10⁻⁵ mM (milli-molar, 1 molar corresponding to a solution of 1 mol/L of concentration) which is a mean oceanic value (Lamontagne et al., 1973).

* Condition 4: sedimentation-related overpressure calculated in one dimension with SeCoV3 is imposed at the southern border of the model, where the sedimentary column is the thickest.

* Condition 5: dissolved-methane concentration imposed at the base of the fault network with values ranging from 56 to 5550 mM, methane being generated in the Maykop Formation, the regional source rock in the SCB.

Several tests were run with different hydraulic diffusivities ($D_h$) of the fault network, considering faults as either transmissive or sealing for lateral fluid flow based on regional studies and previous models (Battani et al., 2010; Bredehoeft et al., 1988; Caine et al., 1996; Evans et al., 1997; Gautherot et al., 2015; Gordon & Flemings, 1998; Javanshir et al., 2015; Wibberley et al., 2017).

### 2.2. Mud ascent

In order to model the mud flowing towards the surface, a fluid mechanics approach is commonly adopted using simplified Navier-Stokes equations (Collignon et al., 2018a; Gisler, 2009; Zoporowski & Miller, 2009). In the present work, we combined two approaches to model the ascent of low-density mud: one solving the problem in two dimensions at low viscosities (Tryggvason, 2012) and one for extrapolating the results to realistic higher viscosities using one-dimensional calculations (Furbish, 1997).

#### 2.2.1. Equations

2D simulations are based on Navier-Stokes equations governed by two basic equations. First, neglecting surface tension, taking gravity as sole body force and assuming constant and uniform viscosity, the momentum equation is given by Tryggvason et al. (2006) as:
$$216 \qquad \rho \frac{\partial u}{\partial t} + \rho \nabla \cdot u u = -\nabla p + \rho g + \mu_0 \nabla^2 u \qquad (2)$$

where $\rho$ is the density of the mud (kg/m³), $u$ is the fluid velocity (m/s), $g$ is the gravitational acceleration (m/s²), $\mu_0$ is the viscosity considered as constant and uniform (Pa.s) and $p$ is the fluid pressure (Pa) (Tryggvason et al., 2006).

Secondly, mass conservation for incompressible fluids is given by:

$$220 \qquad \nabla \cdot u = 0 \qquad (3)$$

as stated in Tryggvason et al. (2006). A Boussinesq approximation is applied to the Navier-Stokes equations to account for mud density variations during its ascent (Boussinesq, 1897). Gas solubility (Henry's law; Henry, 1803) and compressibility (Boyle's law; Boyle, 1965) are used to derive at each time step the effect of gas saturation on the mud density through the following equation:

$$225 \qquad \rho = \phi(1 - S_g)\rho_w + (1 - \phi)\rho_s \qquad (4)$$

where $\phi$ is the porosity and $\rho$, $\rho_w$ and $\rho_s$ are the respective mass-densities of mud, water and solid particles. The density of free gas is considered as being negligible compared to the water density and the density of the solid particles.

Finally, extrapolation to more realistic viscosity values was completed by a simple 1D calculation following a simple example of Furbish (1997) modeling a buoyant magma flow within a vertical dyke. It was possible to calculate the velocity of the ascending mud using the following equation:

$$235 \qquad v(x) = \frac{1}{2\mu}(\rho_c - \rho)g(r^2 - x^2) \qquad (5)$$

with $v$ the velocity of the ascending mud (m/s), $\mu$ the viscosity of the mud (Pa.s), $\rho$ and $\rho_c$ being the respective mass-densities of the mud and of the surrounding rock (reference) (kg/m³), $g$ the acceleration of gravity (m/s²), $r$ the radius of the mud chamber (m) and $x$ the position relative to the center of the mud chamber (m). The density of the reference rock was calculated just above the mud generation zone using the porosity obtained through the 2D sedimentation models and a density for solid grains of 2650 kg/m³. The degree of gas saturation was taken at 38 %, i.e. the limit for mud generation based on the laboratory testing by Blouin et al., (2019b). It allowed the calculation of the mud density using Eq. (4). Therefore, applying Eq. (5), a maximum mud velocity was calculated for different sediment viscosities.

### 2.2.2. Numerical resolution

For 2D simulations, we used Tryggvason et al. (2006) resolution method of the Navier-Stokes equations (Eq. (2) and (3)) by including the effect of free gas expansion on the mud density (Eq. (4)). More precisely, the software used in this paper is based on code1 of Tryggvason (2012). The code is freely available and free of use (Tryggvason, 2011).

The momentum equation is integrated with time by splitting Eq. (2) into a velocity term and a pressure term, integrated separately. Then, each term of Eq. (2) is discretized using a Finite-Volume approach where both equations are applied over a small control volume. More details concerning the basic Navier-Stokes equations and the numerical methods are provided in Tryggvason (2012).
The degree of gas saturation and the and the mud density were calculated in parallel at each calculation time step and introduced as an input to the Navier-Stokes solver. That leads to density inversion between the mud generation zone and the host rock (Brown, 1990; Collignon et al., 2018a; Deville, 2009; Kopf, 2002). Mud and host rocks are modelled as having the same viscosity that does not depend on the degree of gas saturation.

### 2.2.3. Geometry and limit conditions

The 2D simulation of mud ascent was achieved over structural models with 100 x 100 nodes representing a 10-km-long and 4-km-thick section, corresponding to the approximate length of the Absheron anticline in a SSW-NNE direction and to the depth below sea level of the mud generation zone. The water column is considered and serves as an upper boundary limit where no exchange with the mud and sediment zones is allowed. The two lateral boundaries are considered as no-flow boundaries. The initial density of the sediment decreases upward from 2100 kg/m³ at the base of the model to 1900 kg/m³ at the seabed. An initial mud body is considered at the initial state and serves as a constant mud source in the model. It corresponds to the volume of damaged sediments calculated from the regional diffusion simulation and has the same density as unconsolidated marine sediments (1900 kg/m³). This initial volume of mud can be also considered with an initial overpressure equivalent to the one obtained from diffusion simulations.

The mud generation process involves the formation of hydrofractures (Roberts et al., 2010). These hydrofractures should propagate vertically within the extrados of the Absheron anticline where normal faulting was observed. Therefore, vertical conduits having the same density as the mud source can be introduced in the initial geometry in order to simulate the potential presence of vertical fractures.

Viscosities measured for surface marine sediments highly depends on the type of sediments, but are clearly below $10^8$ Pa.s (Jeong, 2013; Locat & Demers, 2008; Torrance, 2010). However, this range of values is not representative of highly consolidated and stratified sediments, whose typical values of viscosity exceed $10^{17}$ Pa.s, value corresponding to evaporitic sediments (Collignon et al., 2018a; Mukherjee et al., 2010). With the depth range and type of host sediment encountered here, taking realistic values would have led to exceedingly long computation times.

## 3. Results

### 3.1. Mud generation

#### 3.1.1. One-dimensional sedimentation and pore pressure accumulation

The Anhydritic Surakhany Formation (ASF) acts as an efficient seal, as overpressure rapidly builds up below this interval at well locations (Blouin et al., 2019a). The layer is assumed to have constant porosity and with a constant low permeability three orders of magnitude lower than the most permeable layer (Layer 4, Figure 2). Its sedimentation rate was taken equal to the one used for Layer 4, as being part of the same overall stratigraphic interval (Figure 2).

Our in-situ pore pressure data show that the presence of the low permeability ASF has prevented fluids from being expelled during the burial of the Productive Series (PS). The increase in the overpressure gradient below the ASF clearly shows that fluids are trapped below the ASF, building-up the overpressure in the underlying strata. Therefore, the Upper PS might have kept an
abnormally high porosity through compaction disequilibrium (Osborne & Swarbrick, 1997). Pressure data reported in Javanshir et al. (2015) as well as the reservoir overpressure values detailed in Blouin et al. (2019a) show that the Middle Productive Series and Lower Productive Series are drained through laterally continuous and connected reservoirs, which should have limited under-compaction in these intervals. Thus, in the following simulations Layer 4 permeability was considered higher than in underlying and overlying intervals using a modified permeability law (Eq. (6)). That modified law reproduces the under-compaction caused by the ASF as displayed in Figure 3:

$$ \ln(K) = 3.06 e - 17.66 \qquad (6) $$

with K the hydraulic conductivity in m/s, and e the void ratio.

Using permeabilities obtained from oedometer tests and well data analysis (Blouin et al., 2019a), initial 1D properties were obtained from the deeper part of the section presented in Figure 2. The considered sedimentation rates are summarized in Figure 3.

Figure 3 displays the main results of the 1D modeling that are taken as boundary condition 4 (Figure 2) in the 2D-diffusion models presented below. The results show the impact of the low-permeability ASF on overpressure that rises sharply when crossing this interval, accounting for the observed pressure build-up in well pressure logs (Blouin et al., 2019a). The low-permeability zone is well expressed in Figure 3 and contrasts with the high-permeability Layer 4. The overpressure slowly builds up over Layer 1, stays constant over Layer 4 and slightly increases over Layer 3, 2 and 1. The highest overpressure reached 40 MPa, which is in accordance with the pore pressure measured at well locations.

### 3.1.2. Two-dimensional transient-diffusion processes, gas exsolution and damage

Four numerical simulations were conducted with different input parameters. The parameters used during the simulations are synthetized in Table 1. Most of the parameters are common for all the simulations and were obtained from literature and industrial reports (seafloor temperature, Henry methane constant) or from laboratory testing and well log interpretation (temperature gradient, sediment petrophysical properties; Blouin et al., 2019a). Four parameters were modified in order to fit observations and regional geology: fault hydraulic diffusivity, fracture condition, starting date for overpressure transmission, initial methane concentration.

The goal of this 2D-modeling is to simulate the mud generation process through sediment damage caused by gas exsolution at the Absheron mud volcano location. Depth is given as meters below mean sea level in this section.

#### i. Case 1: integration of the ASF layer properties

The modeled geometry is displayed in Figure 2. The fracture condition $\Delta u/\sigma'_v$ was initially set at 0.9 given the mean fracture condition obtained from well data (Blouin et al., 2019a). Lateral overpressure transmission started 1 My after inception of methane diffusion in fault networks. Results of this simulation after 5 My are displayed in Figure 4. Methane diffusivity was set at $1.49 \times 10^{-7} \text{ m}^2/\text{s}$ in sediments, and faults have a hydraulic diffusivity of $5.5 \times 10^{-6} \text{ m}^2/\text{s}$ and a relatively high methane diffusivity coefficient of $5 \times 10^{-4} \text{ m}^2/\text{s}$ in order to artificially include gas advection impact and create a preferential advective pathway for fluids.

Figure 4a shows overpressure diffusion across the structural model. Overpressure strongly
increases below the ASF, reaching 38 MPa in the south of the model, and stays lower than 18 MPa above the ASF. Moreover, as Layer 4 is more permeable than the others, overpressure propagates more rapidly along this stratigraphic layer. Figure 4c shows the overpressure build-up below the ASF at the Absheron location that reaches 35 MPa over less than 1000 m.

However, in terms of $\Delta u/\sigma'_{v}$ distribution, Figure 4b clearly shows that the maximum is reached just below the ASF, in the northern edge of the model where $\sigma'_{v}$ is low due to thin overburden and to the boundary condition (Figure 2), which prevents lateral exchanges between the model and the outsides. Hydrofracturing conditions are reached as $\Delta u/\sigma'_{v}$ is above 0.9 (green area in Figure 4b). No fracture conditions are reached at the Absheron location where calculated values do not exceed 0.7. Methane circulation is effective around the fault network showing that over 5 My, sediments get saturated with dissolved methane through the simulated transmissive fault network.

Thus, the presence of the low permeability ASF allows reproducing pressure transmission and build-up mainly below this interval as well as reaching hydrofracturing conditions along the bottom end of the ASF. However, fractures occur at the northern end of the structural model where $\sigma'_{v}$ is lowest and where overpressure accumulates due to model boundary effects. Besides, no free gas, and consequently no sediment damage was observed due to fracture formation in an area where no dissolved gas was present.

## ii. Case 2: sealing faults

Faults were initially set to create preferential pathways for fluids. Nevertheless, faults could be vertical pathways but horizontal seals for fluid circulation due to permeability anisotropy (Caine et al., 1996; Evans et al., 1997; Wibberley et al., 2008). Indeed, different authors noted in several regions around the world that faults may be sealing in the perpendicular direction but transmissive along the fault direction and that their permeability may vary depending on external parameters such as pore pressure, deformation style and the lithologies put in contact through the fault and within the fault itself (Henry et al., 2019; Morley et al., 2017; Wibberley et al., 2017). This feature was also already highlighted over the SCB and particularly in the Absheron anticline where gas was found in the northern flank while no economic gas accumulation was encountered in the southern flank (Gautherot et al., 2015; Grosjean et al., 2009; Javanshir et al., 2015). Several regional fluid flow numerical models already integrated fault anisotropy (Deville et al., 2010; Gordon & Flemings, 1998; Schneider et al., 2004).

To simulate the horizontal sealing effect of the fault, the horizontal hydraulic diffusivity of the fault was taken equal to $5.5 \times 10^{-8} \text{ m}^2/\text{s}$. Vertical transmissivity was simulated by keeping a methane diffusivity of $5 \times 10^{-4} \text{ m}^2/\text{s}$. Results of this simulation are displayed in Figure 5.

Figure 5a displays the 2D distribution of overpressure after 5 My of calculation. The main feature is that with the sealing faults, overpressure builds up south of the faults, reaching 40 MPa, while in the northern part, the overpressure stays limited (around 18 MPa). The effect of the high permeability layer 4 contrasting with the low permeability ASF is still displayed with overpressure transmitting more rapidly along Layer 4 than in overlying intervals and with limited overpressure buildup in Layer 1. Figure 5c shows the overpressure profile with depth at the Absheron location. The sharp overpressure increase across the ASF and the weak overpressure gradient along Layer 1 are well captured by the model. Besides, an overpressure contrast is visible between Layer 4 and Layer 3 due to the effect of sealing faults, creating a trend similar to that observed in pressure logs
at wells, with pressure horns and peaks (Blouin et al., 2019a).

High values of $\Delta u/\sigma'_{v}$ are now distributed in two main areas. To the north of the structural model where $\sigma'_{v}$ is low and where some overpressure accumulated below the ASF, but most of all just south of the sealed fault network and along the ASF at the crest and along the southern flank of the Absheron anticline (Figure 5b). However, the fracture condition 0.9 is not reached, thus no fractures were created. The highest $\Delta u/\sigma'_{v}$ values reached are comprised between 0.70 and 0.75 (Figure 5b). Methane diffusion is similar to what was observed in the previous simulations as methane diffusivity of faults and sediments were not modified (Figure 5b). Methane concentrations are also the same.

Therefore, the fault network now acts as a seal for lateral pressure transmission and as a vertical pathway for fluids, allowing to approach critical conditions at the Absheron location.

### iii. Case 3: new assessment of fracture conditions

The fracture condition was previously set at 0.9 because of the mean value obtained on well data of the ratio between fracture pressure and overburden pressure. However, this value may be as low as 0.8 depending on the considered depth. Besides, at the time when the AMV started being active, sediments may have been less consolidated than at present which led to lower tensile strength hence lower hydrofracturing conditions. Therefore, a fracture condition of 0.7 was considered as plausible and is now considered in the following calculations. This value $\Delta u/\sigma'_{v}$ was already given as a minimum value for brittle failure of sedimentary rocks in several studies (Grauls & Baleix, 1994; Sibson, 2003).

A new calculation only changing the critical fracture criterion from 0.9 to 0.7 was carried out. Figure 6 synthetizes the results after 5 My of calculation.

Figure 6a and Figure 6c are roughly the same as Figure 5, as no changes were made that could affect directly the pressure field. Figure 6b shows that fracture condition was reached and that a fracture opened in the area of the Absheron anticline, south of the fault network, where $\Delta u/\sigma'_{v}$ is the highest. As fractures have a methane diffusivity of $5 \times 10^{-4} \text{ m}^2/\text{s}$, pore fluids migrated through the fractures along a horizontal plane, changing the distribution of the concentration of dissolved methane (Figure 6b). However, no gas exsolution occurred in this model, since the degree of gas saturation stayed at zero. This could be either due to the fact that the pressure drop through fracture was not high enough, or that the dissolved methane concentration is not sufficient, two parameters that control gas exsolution processes (Brown, 1990; Duan & Mao, 2006). Figure 6a and Figure 6c do not display any pressure decrease at fracture depth. As fractures are generated in an area where overpressure is constant, no drainage of fluids was possible, preventing pore pressure from a significant decrease. A vertical fracture that would have connected Layer 4 and Layer 6 through the ASF would have permitted a strong pressure decrease. Nevertheless, this was not possible since the ASF layer is characterized by a low-permeability coefficient impeding vertical pore pressure diffusion and therefore hydrofracturing.

Another important point is that hydrofracturing happened only 500,000 years after the beginning of pressure diffusion calculation, so 1.5 My after the beginning of the simulation. This shows that there is no need for calculating diffusion for a period as long as 5 My and that diffusion calculation could start at 3 Ma, date when all the PS were already deposited (Abreu & Nummedal, 2007; Forte & Cowgill, 2013; Green et al., 2009; Morton et al., 2003; Vincent et al., 2010). This could account for the deposition time of the main reservoirs, when hydrocarbon generation was
already active (Guliyev et al., 2011; Inan et al., 2002; Smith-Rouch, 2006; Tagiyev et al., 1997).

As a conclusion of this calculation, fractures were generated at the Absheron anticline crest, south of the fault network, and dissolved gas diffused into the fracture network. However, no free gas was formed because of either unsubstantial pressure decrease or low dissolved methane concentration.

#### iv. Case 4: high gas concentration at the base of the fault network

In order to test the impact of free gas at the fracture location (Figure 6b), the dissolved methane concentration of the fluid source at the bottom of the fault network (coming from a deeper source such as the Maykop source rock) was increased to 555 mM and to 5550 mM. The latter high methane concentration (twice higher than methane solubility at the entry point temperature and pressure conditions) implies the presence of free gas in addition to dissolved methane. However, this initial free gas volume accumulates during several thousands of years along with sedimentation and consolidation of host sediments though overburden. Thus, it does not cause any damage to the sediment as its volume decreases and as it gets dissolved as overburden increases. Here, damaging is considered to occur only during short-term exsolution process. The simulations start calculating pressure diffusion after 3 My of methane diffusion through the faults and the stratified sediment layers.

The simulation considering an initial methane concentration of 555 mM was not conclusive as the results obtained are the same in Figure 6. However, this showed that pore-pressure diffusion starting at 3 My does not impact the previous results.

Figure 7 shows the results of the simulation with an initial methane concentration of 5550 mM. Figure 7a shows that the overpressure field does not change compared to previous simulations. Figure 7b shows a similar $\Delta u/\sigma'_v$ distribution than in Figure 7b with a fracture generation below the ASF, south of the fault networks. The dissolved methane concentration follows faults and fractures but with higher concentrations, the lower value plotted being 500 mM (Figure 7b). Figure 7c displays the degree of free gas saturation ($S_g$) calculated across the structural model. $S_g$ is above zero in an area corresponding to the fractured area close to the fault network and just below the ASF limit. Maximum values of 1 are reached in the center of the fractures. The model triggers gas exsolution locally after fracture generation for sufficient dissolved methane concentration.

Gas exsolution has a twofold effect. First, dissolved methane distribution is disturbed and partly depleted around the faults compared to the simulation proposed in Figure 6 (see Figure 7b) due to the fact that, locally, a part of the initially dissolved methane exsolved to form free gas. Moreover, as the simulation calculates the preconsolidation pressure ($\sigma'_p$), based on Eq. (1) (Blouin et al., 2019b) linking $S_g$ with a preconsolidation ratio, we observe a local decrease in preconsolidation pressure in response to the presence of free methane (Figure 7d). This decrease in $\sigma'_p$ is observed in the area where $S_g$ is non-zero. $\sigma'_p$ reaches almost 0 kPa where $S_g$ is greater than 0.38, in areas where it was around 50 MPa before gas exsolution (Figure 7d).

Thus, this simulation considering 1- the presence of the ASF; 2- sealing faults; 3- a fracture condition of 0.7; and 4- an initial methane concentration of 5550 mM, was able to model sediment damage in an area just below the ASF and at the Absheron anticline crest, through fracture opening and the subsequent gas exsolution.
The integration of a low permeability layer corresponding to the ASF allowed transmitting overpressure into the Upper Productive Series and reaching critical overpressure conditions. The modification of fault hydraulic diffusivity in order to create horizontal seals resulted in critical conditions at the crest of the Absheron Anticline. These critical conditions created fractures that, coupled with sufficiently high dissolved methane concentration, allowed gas exsolution. Gas exsolution locally decreased the preconsolidation pressure of sediments sufficiently for them to lose all their initial structure. The observations made during the laboratory testing led by Blouin et al. (2019b) allow concluding that with the parameters used, conditions for mud generation are reached in the Upper Productive Series, at the Absheron mud volcano location.

## 3.2. Mud ascent

The approximate volume of mud and the initial geometrical conditions for mud ascent modeling were estimated from the previous calculations of the mud generation (Figure 7). The initial (2D) mud body for all simulations was set as a convex-up half-disc with a radius of 500 m, its base (diameter of the half-disc) being located 3400 m below seafloor with 500 m of water column above and centered on the horizontal axis of the model (Table 2). When considered, the maximum initial mud overpressure was set at 35 MPa, which is the pressure obtained in the mud generation zone at the end of the simulation presented in Figure 7. Mass-densities used for all simulations are indicated in Table 2.

Three values of the common viscosity of the ascending mud and host sediments were tested as well as the impact of the presence of vertical conduits and of their length (Table 2).

### 3.2.1. 2D simulation results

Figure 8 gives an example of the 2D simulation results in terms of density distribution, sediment velocities and degree of gas saturation. This particular example displays the case where $\mu = 10^9$ Pa.s without an initial vertical conduit and no overpressured mud. From this example the main steps of each simulation can be described:

1. Initial geometry with the source mud chamber having a density of 1900 kg/m³ and with the rest of the sedimentary column ranging from 2100 kg/m³ at 3900 m to 1900 kg/m³ at the seafloor. Initial degree of gas saturation is 0.38, the value for which mud generation happened in 2D mud generation simulations (Figure 7 and section 3.1.2.iv). Sediment velocities are equal to zero at first.

2. Initiation of vertical mud ascent after an initial stability. Sediments above the ascending column starts having positive velocity values. The density remains between 1900 and 1950 kg/m³ and the gas saturation increases slightly towards the top of the ascending column.

3. Acceleration of the mud ascent as the degree of gas saturation increases in the ascending mud mass. Increase in sediment velocities and in the length of the affected strata.

4. The ascension slows down as the ascending mud gets close to sediment/water interface acting as a no-exchange boundary. The ascending mud top flattens. This is considered as the final state where extrusion would occur. Sediment velocities reach their maximum values along the ascending mud column.
Results from the different simulation cases (Table 2) are compared in Table 3 in terms of upward movement initiation time, extrusion time, maximum sediment velocities and final diameter of the ascending mud column. These results highly depend on the spatial and temporal resolution of the simulations, hence small variations are not considered to be relevant. Mud ascent in simulations with higher viscosities (e.g. $\mu = 10^{10}$ Pa.s; Table 2) never initiated after 24 hours of calculation. Simulation results shown in Table 3 indicate that when viscosity increases by one order of magnitude, velocity decreases by one order of magnitude while the timing of initiation and extrusion increases by one order of magnitude (Table 3). When mud overpressure is considered in the simulations, it decreases the time for initiation and extrusion by a ratio of 0.65 to 0.5 (Table 3). When mud overpressure is not considered, the presence of vertical conduits above the mud chamber mainly influences the $V_{max}$ which seems to decrease as the length of the conduit increases (Table 3). For overpressured muds, $V_{max}$ seems to remain roughly constant. The presence of a conduit also influences the final diameter of the gassy-mud column that decreases as the length of the conduit increases (Table 3).

### 3.2.2. Extrapolation to realistic viscosities: 1D calculations

In order to extrapolate the 2D simulation results to higher viscosity values, we tested simple 1D calculations based on the case of a buoyant magma flow along a vertical dyke (Furbish, 1997) in order to compare the velocities and timing obtained. Using Eq. (5) with r the radius of the mud chamber ($r = 500$ m) and taking $x = 0$ (maximum velocity above the center of the mud chamber) it was possible to obtain the maximum velocity above the mud chamber for varying viscosities. Figure 9a shows the results of 1D calculations of maximum velocities ($V_{max}$) for viscosities comprised between $10^{-8}$ and $10^{15}$ Pa.s, considered as a valid approximation range for mud viscosities that are typically below the viscosity of sedimentary rocks. Two sets of models were run: the first one without initial mud overpressure (black), corresponding to the mud density calculated through Eq. (4); the second (red line) considers an initial mud overpressure. Overpressure is limited to 78% of the effective stress (value of $\Delta u/\sigma'_v$ near the mud source at the end of 2D diffusion models; Figure 7), thus, the density of mud was multiplied by 0.78 to take the mud overpressure into account. The maximum velocities obtained above the mud chambers for the 2D simulations without conduit are also plotted in Figure 9a. Figure 9b displays the same results in terms of minimum time for extrusion that was calculated considering that the velocity is $V_{max}$ all along the conduit. The results in Figure 9b clearly show that the minimum time needed for the buoyant mud to reach the seafloor and form a mud volcano depends on mud viscosity.

One-dimension calculations show that the maximum velocity is inversely proportional to the viscosity (Figure 9a). For constant viscosity, $V_{max}$ increases by a factor 4.5 when considering an initial mud overpressure in the calculation. Thus, the two trends, with and without mud overpressure, are parallel. For higher viscosities ($10^{15}$ Pa.s), the maximum velocity is comprised between $10^{-6}$ and $10^{-7}$ m/s. Logically, the time for extrusion increases by an order of magnitude for each order of magnitude of $\mu$ (Figure 9b), a result which directly arises from the form of Eq. (5). The consideration of the mud overpressure allows reducing this time for extrusion. Therefore, for $\mu = 10^{15}$ Pa.s, it would take a minimum of 500 years for the buoyant mud to reach the surface if mud overpressure was not considered, while it only takes 100 years if mud overpressure is integrated in the calculation.

The 2D simulations results display the same relationship. $\mu$ and $V_{max}$ are inversely proportional (Figure 9a). Considering the initial mud overpressure leads to a 4-fold increase of
Vmax compared with the case without overpressured mud (Figure 9a). Moreover, the values obtained without overpressured mud stands close to the black line (1D calculations without mud overpressure), and the values considering mud overpressure are closer to the red line obtained from simulations without mud overpressure (1D calculations considering mud overpressure).

The same calculation was applied but in the case of a 50-m-wide conduit corresponding to the size of the conduits in the 2D simulations (Figure 10). The same relationship between μ and Vmax as well as μ and the time for extrusion, are displayed in Figure 10. The impact of mud overpressure seems to be limited as the difference between the Vmax and the time for extrusion obtained considering mud overpressure and without mud overpressure is reduced compared to Figure 9. Moreover, Vmax values extracted from 2D simulations are larger by nearly one order of magnitude than those obtained from 1D calculations (Figure 10a). Therefore, the time for extrusion is smaller by one order of magnitude (Figure 10b). Besides, the conduit width corresponds to the minimum that could be computed considering the model resolution, one node corresponding to 100 m and the conduit being computed over a unique node. Hence, the results displayed are not representative of the circulation that would occur in natural fractures whose width is much smaller.

From this 1D calculation, it arises that the time required for achieving extrusion, considering a viscosity of 10<sup>15</sup> Pa.s, is significantly increased. Considering mud overpressure, it would take more than 10,000 years for the mud to reach the surface through thin fractures when it was 100 years for the calculation considering the mud chamber radius.

## 4. Discussion

### 4.1. Synthesis of the main results

Representing the low-permeability Anhydritic Surakhany Formation (ASF) in the 1D model generates a sedimentation-related overpressure contrast between layer 6 and deeper layers, thus reproducing the sharp overpressure increase observed in Absheron well-log data (Blouin et al., 2019a). Moreover, the model simulates maximum overpressure values of 40 MPa corresponding to the values measured in situ at the Absheron exploration wells. It also improved the preliminary 2D diffusion model presented in Blouin et al. (2019a) with overpressure occurring essentially below the ASF and transmitted through Layer 4 (Figure 4). Overpressure diffusion started after 3 My of simulation, a duration corresponding to the deposition of the entire Productive Series interval. Lateral compartmentalization of the 2D model through the presence of faults having low horizontal permeability (Gautherot et al., 2015; Javanshir et al., 2015) was necessary to generate overpressure at the location of the AMV (Figure 5). Faults are commonly regarded as presenting permeability anisotropy, with the higher permeability along the fault surface (Caine et al., 1996; Deville et al., 2010; Evans et al., 1997; Henry et al., 2019; Morley et al., 2017; Schneider et al., 2004; Wibberley et al., 2008), thus methane-saturated water migration from the Maykop Fm. was provided by high methane diffusivity in faults, artificially simulating high vertical permeability accounting for advection processes. A new fracture condition value of 0.7 based on pressure logs at well sites was adopted, which generated in the model a fractured area at the Absheron crest, just below the ASF (Figure 6). For an initial methane concentration of 5550 mM, gas exsolution was triggered in the fractured area resulting in a strong sediment damage in the fracture zone, with preconsolidation pressure locally approaching zero (Figure 7). From the results of the laboratory testing discussed by Blouin et al. (2019b) the improved 2D diffusion model simulated the mud generation zone formation below the Absheron MV.
2D-simulation of mud remobilization show that ascent velocity is inversely proportional to mud viscosity, while the time needed for mud ascent initiation and for mud extrusion is proportional to viscosity (Table 3). When mud overpressure is considered in the simulations, it decreases the time for initiation and extrusion by a ratio of 0.5-0.65 (Table 3). The presence of conduits mainly influences the final width of the gassy-mud column, which decreases with an increasing conduit length (Table 3). Through one-dimensional extrapolation it would take a minimum of 500 years for a buoyant high-viscosity mud to reach the surface if mud overpressure is not considered, while it would take only 100 years if the mud overpressure is integrated in the calculation (Figure 9b). Considering the radius of a conduit (Figure 10b), it arises that the timing for extrusion is significantly increased as it would take more than 10,000 years for the same overpressured mud to reach the surface. Combining 1D and 2D simulations, adapting the Navier-Stokes equations to consider the impact of gas expansion on sediment density through Boyle's law and the influence of mud overpressure over the ascension velocity allowed estimating the time needed for mud generated at depth to reach the seafloor and to form a mud volcano through density-inversion.

## 4.2. Overpressure, hydrofracturing and gas exsolution: from stratified sediments to fluid mud

Calculations of pressure related to the maximum gas column that could have been trapped within the ASF (chapter 3, Blouin, 2019) showed that the sole presence of a gas reservoir at the crest of the Absheron fold cannot explain the formation of the Absheron mud volcano, as the fracture pressure was not reached. Hydrofracturing is believed to be necessary to initiate a MV formation as demonstrated through the monitoring of the Lusi initiation (Tingay et al., 2017) and through the analysis of the Absheron Mud Volcano dataset (Blouin et al., 2019a). Therefore, additional parameters are required to explain the AMV formation. The sedimentation rates recorded in the SCB are among the highest in foreland basins (Tagiyev et al., 1997; Allen et al., 2002; Smith-Rouch, 2006; Egan et al., 2009; Green et al., 2009). High sedimentation generates overpressure in low-permeability strata (Dugan & Flemings, 2000; Opara, 2011; Osborne & Swarbrick, 1997) and overpressure is one of the key triggers and drivers for MV formation (Deville, 2009; Dimitrov, 2002; Kopf, 2002; Mazzini & Etiope, 2017).

The 2-D diffusion model presented in this paper is based on regional data extracted from literature (Diaconescu et al., 2001; Green et al., 2009), on observations and in situ measurements made on the Absheron gas condensate field, as well as on the laboratory experiments conducted by Blouin et al. (2019b). The regional-scale pressure gradient and northward fluid circulation from the deep SCB was already observed and discussed by several studies (Bredehoeft et al., 1988; Gautherot et al., 2015; Grosjean et al., 2009; Javanshir et al., 2015), and this type of modeling, studying the influence of regional physical parameters such as sedimentation rate and pressure gradient on local structures such as landslides or venting sites, has already been applied to other sedimentary basins (Dugan & Flemings, 2000; Schneider et al., 2004; Kvalstad et al., 2005; Hustoft et al., 2009). Deville et al. (2010) applied a similar fluid flow numerical model to show that MVs in Trinidad are located above overpressured areas. However, it is the first time, in our knowledge, that a 2D-diffusion model is applied to simulate sediment remobilization into mud at depth. The 2D modeling presented in this study allows drawing a coupled process for mud generation below the AMV. Due to high sedimentation rates, overpressure is generated in the deeper part of the SCB and is transmitted northwards along under-compacted layers that kept an abnormally high permeability. The fault system at the core of the Absheron anticline leads to local
overpressure build-up above the crest of the anticline, eventually initiating hydrofracturing and saturating the pore waters of the Upper Productive Series with methane. Hydrofracturing may trigger a local decrease in overpressure, allowing exsolution of the dissolved methane and triggering sediment damage when S$_{g}$ exceeds 0.38 (Blouin et al., 2019b). Thus, the 2D diffusion model integrating the mechanical properties of sediments, sedimentation rates, compartmentalization effect of faults and regional seal and the impact of gas exsolution reproduces mud generation conditions at depth below the AMV.

### 4.3. Mud extrusion resulting from density-inversion driven by gas exsolution

Once mud is generated at depth, it is necessary to understand how it is transferred to the surface. Modeling shows that the mud ascent time is proportional to the mud viscosity. Taking mud overpressure into account improves the convergence of 1D and 2D results (Figure 9). Similar values of velocity were obtained by Collignon et al. (2018b) who calculated flow rates depending on pipe radii, particle size and density and the type of gas for equivalent viscosities. One-dimensional extrapolation allows to test viscosities closer to those of sedimentary rocks, which was not possible on the 2D simulations due to time limitations. It would take approximately 100 years for the mud generated at the AMV location to rise to the surface if mud was only transported in response to density inversion and mud overpressure. The link between 2D results and the extrapolation in 1D may not be straightforward. The extrapolation in 1D is made over 6 orders of magnitude of viscosities and compared to only two 2D control points. The relationship between mud remobilization and the host sediment viscosity may not be a linear process when dealing with high values of µ. Therefore, the results have to be addressed with caution. Nevertheless, some MVs display very slow extrusion rates, such as the Kotyrdag (onshore Azerbaijan), that extrudes high viscosity and slightly wet mud with an estimated rate of 0.02 m/day (Dupuis, 2017), in other words 2.3 x 10⁻⁷ m/s, which is the order of magnitude obtained from our modeling.

The sole effect of gas expansion on the mud density is thus able to drive mud ascent up to the seafloor. The process may be accelerated if the influence of gas exsolution and expansion on mud viscosity and of fracture propagation ahead of the rising mud are accounted for. This modeling work resulted in extrusion rates of the same order of magnitude that what has been observed on MVs extruding high viscosity mud.

### 4.4. Towards a quantitative formation model for the Absheron mud volcano

The conceptual formation sketch for the AMV presented in Blouin et al. (2019a) was improved by integrating the main findings of Blouin et al. (2019b) and of the numerical models carried-out in the present study. The updated formation sketch is composed of 7 different phases shown in the different panels of Figure 11:

Phase 1 - At 3.4 My, end of the rapid deposition of the Productive Series (over 3.5 km in 2 My) above the gas-mature Maykop Formation that generated hydrocarbons diffusing slowly through the sedimentary column (Figure 11a). The Maykop Formation started generating gas during the Late Miocene in the Shah Deniz region located in a structural setting comparable to Absheron (figure 2.9 in Alizadeh et al., 2017). Moreover, the deposition of the Anhydritic Surakhany Formation provided an efficient seal allowing a slow and uniform overpressure buildup in the Productive Series generated by high sedimentation rates.

Phase 2 - At 0.7 Ma, the growth rate of the Absheron fold increased significantly. Methane
circulation was therefore focused into the faulted core of the anticline (Figure 11b). The difference of overburden between the anticline and adjacent syncline generated an overpressure gradient leading to overpressure build-up against the sealing faults (Figure 11b). Besides, at the anticline crest, the vertical effective stress $\sigma'_{v}$ is lower than on the flanks due to reduced overburden thickness. Thus, overpressure generation and transmission as well as low $\sigma'_{v}$ at the crest increased the $\Delta u/\sigma'_{v}$ ratio.

Phase 3 – At $\approx 0.6$ Ma, $\Delta u/\sigma'_{v}$ reached the hydro-fracturing threshold of 0.7 (Figure 11c), allowing hydrofracturing. This led to a local decrease of overpressure through fluid pressure dissipation and triggered methane exsolution where the dissolved-methane concentration reached the methane solubility (Henry's law; Figure 11c). For a gas saturation higher than 0.38, the sediment lost its preconsolidation pressure leading to mud generation (Figure 11c).

Phase 4 – The presence of gas bubbles in the mud significantly reduced its density and its viscosity (Figure 11g). Therefore, the gassy mud started to rise along the open fractures. As gas bubbles moved up, they expanded (following Boyle's law), hence increasing gas saturation and decreasing the mud density and viscosity (Figure 11g). Thus, mud velocity increased as bubbles expanded. Gas expansion maintained the overpressure as mud rose along fractures, leading to a critical value of $\Delta u/\sigma'_{v}$ and to fracture propagation towards the surface (Figure 11g). This stage should have happened over a short period of time, ranging from hours to days and is based on field observations from Roberts et al. (2010) who show that feeding pipes are composed of intensely fractured host rock.

Phase 5 – After several days, the fractured “pipe” connected the mud generation zone, where gas exsolution was still damaging sediments, with the seafloor (Figure 11d). First mud extrusion happened at the seafloor leading to the progressive formation of the extrusive edifice, while mud slowly degassed in the water column (Figure 11d). This direct hydraulic connection prevented overpressure from building up again at depth as it was directly transmitted to the surface (Figure 11d), thus preventing gas exsolution from stopping. At depth, the extrusion of remobilized mud slowly led to roof collapse, progressively forming what is known as a “depletion zone” (Dupuis et al., 2019; Kirkham et al., 2017; Stewart & Davies, 2006).

Phase 6 – After several years, multiple mud extrusion episodes reached the surface as the depletion zone became bigger along with the collapse of the overlying strata, forming a giant mud shield at the seafloor (Figure 11e). The AMV then entered a quiescent phase. Overpressure was not maintained high enough to keep fractures open and overpressure built up again along the sealing faults, preventing further gas exsolution (Figure 11e).

Phase 7 – Present day geometry, after multiple phases of quiescence and activity, creating a complex interdigitated geometry where normal sedimentation predominates during quiescent phases of the volcano. When $\Delta u/\sigma'_{v}$ reaches the fracture opening threshold, further depletion happens in the source, enhancing collapse, while more mud is extruded at the seafloor during the active phases (Figure 11f).

## 4.5. Limitations and potential for improvement

This study and discussions allowed highlighting several limitations and improvements for the two models presented in this paper. For the mud generation simulation:

* More precise measurements of layer permeability and fracture pressure through
721 laboratory testing on rock fragments. This would significantly improve hydro-mechanical properties input for each modeled layer. However, rock cores were only sampled from targeted reservoirs at well sites.

722
723
* Fracturing process could be improved. Fracturing pressure was only estimated from pressure logs (Blouin et al., 2019a). However, LOT/FIT (leak-off tests and formation integrity tests) data indicate that hydrofracturing may occur for an overpressure of 25 MPa (and not at 35 MPa as in Figure 7). That would significantly reduce the simulation time needed to trigger hydrofractures. Besides, vertical fractures generated in the simulation propagates horizontally where a vertical propagation of fractures would be expected. The presence of the low-permeability ASF prevented pore pressure and therefore fractures to propagate vertically even if the maximum principal stress is vertical, the simulation ignoring the impact of compressive tectonics that are observed over the basin.

734
735
736
737
738
739
* The processes being modeled here are 3D processes. Therefore, the results of the simulations are only semi-quantitative. To reach quantitative results, 3D basin modeling should be considered to model fluid migrations accurately (Bruns et al., 2013; Schneider et al., 2000). Faults may act as baffles but are finite features. The geometry of the Absheron anticline being a local topographic high in 3D is also an important feature that 2D simulations does not highlight entirely.

740 For the mud ascent simulation, gas was only regarded as a mean of reducing mud density.
741 However, it could have further impact on mud ascent that are not currently considered by this
742 model:

743
744
745
746
747
748
749
750
* The influence of gas exsolution on mud rheology was not accounted for. The presence of gas and liquids in the mud are known to strongly affect its viscosity (Kopf, 2002; Mazzini & Etiope, 2017). Collignon et al. (2018b) and Zoporowski & Miller (2009) both have considered a dynamic viscosity in their simulations and show that the increasing contrast between mud viscosity and the viscosity of the surrounding rock increases flow rates. Thus, the integration of a dynamic viscosity depending on gas content would certainly increase the velocity of the extrusion process in our model.

751
752
753
754
755
756
757
758
759
760
761
762
* One-dimensional calculations considering fracture radius show that the extrusion time is close to 10,000 years for μ = 10<sup>15</sup> Pa.s. The simulations shown here do not address multi-phase flow, gas being only considered in the density calculation. If free-gas was able to flow directly through the fracture, it may be able to fracture and/or damage shallower sediments, forming a secondary mud generation zone. This was observed at Lusi, where the BJP-1 well transferred overpressured fluids charged with gas to shallower and more fragile strata, triggering hydrofractures and mud generation (Tingay et al., 2008; Tingay et al., 2015). The impact of gas expansion on overpressure and on rock failure is not explored but is expected to reduce significantly the time for extrusion, the timing observed at Lusi being 24-48 hours separating free-gas detection downhole and mud extrusion at the surface (M. Tingay et al., 2017).

763 Further considerations could be assessed in order to improve the mud ascent model:

764
* Field work reported by Roberts et al. (2010) seems to indicate that mud volcano
conduits consist of an intricate and complex fractured volume with detached and rotated blocks. Roberts et al. (2010) explain their observations by the rising of overpressured and low-density mud that generates and propagates fractures as it rises, slowly opening its way up to the surface. Thus, the 2D simulations should rather compute the progressive opening of fractures as mud overpressure increases due to gas expansion. This would imply a combination of solid-rock mechanics and fluid mechanics that represents a complete field of research in modeling and computation (Mainguy & Longuemare, 2002; Wu & Olson, 2015).

* The code of Tryggvason (2012) was developed to solve short-term fluid mechanics problems. Working with values corresponding to a geological background leads to very long calculation time due to computing limitations arising from the resolution method: time increment depends on the mesh resolution and on the viscosity. For instance, for a viscosity of $10^{10}$ Pa.s and for the used mesh of 100 x 100 nodes (mesh resolution of 100 m horizontally and 39 m vertically), a time increment dt of $10^{-5}$ s was needed in order to reach numerical convergence. This implies $1.7 \times 10^{9}$ calculations to simulate the 24 hours of mud ascent. Therefore, another resolution method should be explored, such as implicit time integration that would allow larger time steps (Bijl et al., 2001).

## 5. Conclusion

This study explored the possibility of simulating mud generation and extrusion processes through numerical modeling based on simplified working hypothesis. The ultimate goal is to simulate the complete formation of the Absheron Mud Volcano (AMV) and to quantify the overpressure conditions and methane concentration that led to the present structure. A first model, considering 2D-diffusion laws (Darcy's law and Fick's law) was used to explain the location of the AMV as well as the conditions required to generate fluid mud from stratified solid sediments. The second model applies a fluid mechanics approach to test whether the sole impact of gas expansion on the mud density can lead to mud extrusion. The main results are:

1. The 2D-diffusion model resulted in critical fracture conditions at the crest of the Absheron anticline, below the ASF, where the depleted area was observed by Blouin et al. (2019a). Sedimentation rates to calculate overpressure in the deep basin, considering the low-permeability Anhydritic Surakhany Formation and sealing faults as flow baffles are the key features driving the mud generation process.

2. A methane concentration of 5550 mM, imposed at the base of the fault network considered as permeable in the vertical direction, creates the conditions for triggering gas exsolution and subsequent loss in preconsolidation pressure showing that liquid mud is generated at depth.

3. Mud ascent up to the seafloor through the sole density-inversion provoked by gas expansion is possible over a period of 100 years minimum. This period is in accordance with extrusion rates measured at particular MVs.

4. To accelerate the process, additional parameters such as non-constant viscosity and the upward propagation of fracture as consequences of gas expansion must be considered.

5. The results allow going from a purely conceptual formation model to a semi-quantitative model considering dynamics of the processes involved and quantification of the main controlling physical parameters.
We show that simple physical models integrating realistic geological and hydro-mechanical parameters and behavior of sediments can reproduce the conditions that initiated the formation of the Absheron Mud Volcano. We simulated the key processes inferred from the available dense dataset and from the present geometry described in the literature. These processes are driven by the presence of gas in sedimentary layers coupled with geodynamic compressive regime and excess pore pressure. We demonstrate that a decrease of the formation pore pressure generated by natural (fractures) or anthropogenic (drilling) reasons may cause substantial changes in sediment structure leading to fluid mud generation and the birth of a mud volcano. Therefore, similar modeling work could be undertaken over mud volcanoes located on similar geological contexts where dense geological and geophysical data is available. Such models may help identifying critical conditions of future MVs formation.

## Acknowledgments, Samples, and Data

This study was realized in the scope of the PhD of Arthur Blouin under the joint direction of the University of Pau and Pays Adour (UPPA), Ifremer and Total S.A. To access the SeCoV3 software and output of the different presented simulations and calculations, please directly download it here: https://doi.org/10.6084/m9 figshare.13142993.v1.

## References

826 Abreu, V., & Nummedal, D. (2007). Miocene to Quaternary Sequence Stratigraphy of the South and Central Caspian Basins. *AAPG Studies in Geology*, *55*, 65–86.
828 https://doi.org/10.1306/1205845St553000

829 Alizadeh, A. A., Guliyev, I. S., Kadirov, F. A., & Eppelbaum, L. V. (2017). *Geosciences of Azerbaijan Volume II: Economic Geology and Applied Geophysics* (Vol. II). Basel: Springer
830
831 International Publishing Switzerland. https://doi.org/10.1007/978-3-319-40493-6

832 Allen, M. B., Jones, S., Ismail-Zadeh, A., Simmons, M., & Anderson, L. (2002). Onset of subduction as the cause of rapid Pliocene-Quaternary subsidence in the South Caspian basin. *Geology*, *30*(9), 775–778. https://doi.org/10.1130/0091-7613(2002)030<0775:OOSATC>2.0.CO;2

836 Battani, A., Prinzhofer, A., Deville, E., & Ballentine, C. J. (2010). Trinidad Mud Volcanoes: The Origin of the Gas. *Shale Tectonics: AAPG Memoir*, *93*, 225–238.
837
838 https://doi.org/10.1306/13231317M933427

839 Bijl, H., Carpenter, M. H., & Vasta, V. N. (2001). ATime Integration Schemes for the Unsteady
840 Navier-Stokes Equations. In *15th Computational Fluid Dynamics Conference. Anaheim, CA*.

841 Blouin, A. (2019). *Formation de boue à partir de sédiments stratifiés dans un contexte de volcanisme de boue: le rôle du gaz*. Université de Pau et des Pays de l'Adour.

843 Blouin, A., Imbert, P., Sultan, N., & Callot, J. (2019). Evolution model for the Absheron mud
844 volcano: from in-situ observations to numerical modeling. *Journal of Geophysical
845 Research: Earth Surface, 124*(3), 766–794. https://doi.org/10.1029/2018JF004872

846 Blouin, A., Sultan, N., Callot, J., & Imbert, P. (2019). Sediment damage caused by gas exsolution: a key mechanism for mud volcano formation. *Engineering Geology*, 263.

848 Boussinesq, J. (1897). *Théorie de l'écoulement tourbillonnant et tumultueux des liquides dans les*
849 *lits rectilignes a grande section* (Vol. 1).

850 Boyle, R. (1965). 1662 "A Defence of the Doctrine Touching the Spring and Weight of the Air."

851 *The Works of Robert Boyle, 3, 3–107.*

852 Bredehoeft, J. D., Djevanshir, R. D., & Belitz, K. R. (1988). Lateral fluid flow in a compacting

853 sand-shale sequence: South Caspian basin. *American Association of Petroleum Geologists*

854 *Bulletin*, 72(4), 416–424. https://doi.org/10.1306/703C9A1E-1707-11D7-

855 8645000102C1865D

856 Brown, K. M. (1990). The nature and hydrogeologic significance of mud diapirs and diatremes for

857 accretionary systems. *Journal of Geophysical Research*, 95(B6), 8969.

858 https://doi.org/10.1029/JB095iB06p08969

859 Brunet, M. F., Ershov, A. V, Volozh, Y. A., Korotaev, M. V, Antipov, M. P., & Cadet, J. P. (2007).

860 Precaspian and south Caspian Basins: Subsidence Evolution of Two Superdeep Basins. *Oil*

861 *and Gas of the Greater Caspian Area*, 55(6 mi), 151–155.

862 https://doi.org/10.1306/1205830St553247

863 Bruns, B., di Primio, R., Berner, U., & Littke, R. (2013). Petroleum system evolution in the

864 inverted Lower Saxony Basin, northwest Germany: A 3D basin modeling study. *Geofluids*,

865 13(2), 246–271. https://doi.org/10.1111/gfl.12016

866 Caine, J. S., Evans, J. P., & Forster, C. B. (1996). Fault zone architecture and permeability

867 structure. *Geology*, 24(11), 1025–1028. https://doi.org/10.1130/0091-

868 7613(1996)024<1025:FZAAPS>2.3.CO;2

869 Collignon, M., Mazzini, A., Schmid, D. W., & Lupi, M. (2018). Modelling fluid flow in active

870 clastic piercements: Challenges and approaches. *Marine and Petroleum Geology*, 90, 157–

871 172. https://doi.org/10.1016/j.marpetgeo.2017.09.033

872 Collignon, M., Schmid, D. W., Galerne, C., Lupi, M., & Mazzini, A. (2018). Modelling fluid flow

873 in clastic eruptions: Application to the Lusi mud eruption. *Marine and Petroleum Geology*,

874 90(August 2017), 173–190. https://doi.org/10.1016/j.marpetgeo.2017.08.011

875 Contet, J., & Unterseh, S. (2015). Multiscale site investigation of a giant mud-volcano offshore

876 Azerbaijan - Impact on subsea field development. *Offshore Technology Conference, OTC-*

877 25864(Mv), 1–10. https://doi.org/10.4043/25864-MS

878 DeGroot, D. J., Lunne, T., & Tjelta, T. I. (2010). Recommended best practise for geotechnical site

879 characterisation of cohesive offshore sediments. In S. Gourvenec & D. White (Eds.),

880 *Frontiers in Offshore Geotechnics II* (pp. 33–57). Perth, Australia: CRC Press, Taylor &

881 Francis Group.

882 Deville, E. (2009). Mud Volcano Systems. In N. Lewis & A. Moretti (Eds.), *Volcanoes: *

883 Formation, Eruptions and Modelling* (pp. 95-125 (404)). Nova Science Publishers.

884 Deville, É., Guerlais, S. H., Lallemant, S., & Schneider, F. (2010). Fluid dynamics and subsurface

885 sediment mobilization processes: An overview from Southeast Caribbean. *Basin Research,*

886 22(4), 361–379. https://doi.org/10.1111/j.1365-2117.2010.00474.x

887 Diaconescu, C. C., Kieckhefer, R. M., & Knapp, J. H. (2001). Geophysical evidence for gas

888 hydrates in the deep water of the South Caspian Basin, Azerbaijan. *Marine and Petroleum*

889 *Geology*, 18(2), 209–221. https://doi.org/10.1016/S0264-8172(00)00061-1
890 Dimitrov, L. I. (2002). Mud volcanoes-the most important pathway for degassing deeply buried
sediments. *Earth-Science Reviews*, 59(1-4), 49-76. https://doi.org/10.1016/S0012-
892 8252(02)00069-7

893 Duan, Z., & Mao, S. (2006). A thermodynamic model for calculating methane solubility, density
894 and gas phase composition of methane-bearing aqueous fluids from 273 to 523 K and from 1
895 to 2000 bar. *Geochimica et Cosmochimica Acta*, 70(13), 3369-3386.
896 https://doi.org/10.1016/j.gca.2006.03.018

897 Dugan, B., & Flemings, P. (2000). Overpressure and fluid flow in the new jersey continental slope: implications for slope failure and cold seeps. *Science (New York, N.Y.)*, 289(5477), 288-291.
898 https://doi.org/10.1126/science.289.5477.288
899

900 Dupuis, M. (2017). *Processus de mise en place et évolution des systèmes de volcans de boue*.
901 Université Lille 1 - Sciences et technologies.

902 Dupuis, M., Imbert, P., Odonne, F., & Vendeville, B. (2019). Mud volcanism by repeated roof collapse : 3D architecture and evolution of a mud volcano cluster o ff shore Nigeria. *Marine and Petroleum Geology*, 110(July), 368-387.
905 https://doi.org/10.1016/j.marpetgeo.2019.07.033

906 Egan, S. S., Mosar, J., Brunet, M.-F., & Kangarli, T. (2009). Subsidence and uplift mechanisms
907 within the South Caspian Basin: insights from the onshore and offshore Azerbaijan region.
908 Geological Society, London, Special Publications, 312(1), 219-240.
909 https://doi.org/10.1144/SP312.11

910 Esrig, M. I., & Kirby, R. C. (1977). Implications of gas content for predicting the stability of
911 submarine slopes. *Marine Geotechnology*, 2(1-4), 81-100.
912 https://doi.org/10.1080/10641197709379771

913 Evans, J. P., Forster, C. B., & Goddard, J. V. (1997). Permeability of fault-related rocks, and
914 implications for hydraulic structure of fault zones. *Journal of Structural Geology*, 19(11),
915 1393-1404. https://doi.org/https://doi.org/10.1016/S0191-8141(97)00057-6

916 Forte, A. M., & Cowgill, E. (2013). Late Cenozoic base-level variations of the Caspian Sea: A
917 review of its history and proposed driving mechanisms. *Palaeogeography,
918 Palaeoclimatology, Palaeoecology*, 386, 392-407.
919 https://doi.org/10.1016/j.palaeo.2013.05.035

920 Fowler, S. R., Mildenhall, J., Zalova, S., Riley, G., Elsley, G., Desplanques, A., & Guliyev, F.
921 (2000). Mud volcanoes and structural development on Shah Deniz. *Journal of Petroleum
922 Science and Engineering*, 28(4), 189-206.

923 Furbish, D. J. (1997). *Fluid Physics in Geology*. New York: Oxford University Press.

924 Gautherot, T., Total, S. A., Bakirov, S., Upstream, S., & International, M. (2015). Absheron Field
925 – ABX-2 Exploration Drilling Challenges and Realizations Return to Absheron field, (June
926 2009), 4–6.

927 Gisler, G. (2009). Simulations of the explosive eruption of superheated fluids through deformable
928 media. *Marine and Petroleum Geology*, 26(9), 1888-1895.
929 https://doi.org/10.1016/j.marpetgeo.2008.12.006

930 Golonka, J. (2007). *Geodynamic Evolution of the South Caspian Basin. Oil and Gas of the Greater*
931 *Caspian Area*, **55**, 17–41. https://doi.org/10.1306/1205844St551463

932 Gordon, D. S., & Flemings, P. B. (1998). Generation of overpressure and compaction-driven fluid flow in a Plio-Pleistocene Eugene Island 330, offshore Louisiana. *Basin Research*, **10**, 177–196.

935 Grauls, D. J., & Baleix, J. M. (1994). Role of overpressures and in situ stresses in fault-controlled hydrocarbon migration: a case study. *Marine and Petroleum Geology*, **11**(6), 734–742.

937 Green, T., Abdullayev, N., Hossack, J., Riley, G., & Roberts, A. M. (2009). Sedimentation and subsidence in the South Caspian Basin, Azerbaijan. *Geological Society, London, Special Publications*, **312**(1), 241–260. https://doi.org/10.1144/SP312.12

940 Grosjean, Y., Zaugg, P., Gaulier, J.-M., & Total S.A. (2009). Burial Hydrodynamics and Subtle Hydrocarbon Trap Evaluation: From the Mahakam Delta to the South Caspian Sea. *International Petroleum Technology Conference Held in Doha, Qatar*, 7-9 December 2009, 1–12.

944 Guliyev, I., Aliyeva, E., Huseynov, D., Feyzullayev, A., & Mamedov, P. (2011). Hydrocarbon Potential of Ultra Deep Deposits in the South Caspian Basin. In *AAPG European Region Annual Conference (Vol. 1, p. 66)*. Kiev, Ukraine.

947 Henry, P., Guglielmi, Y., Gout, C., Castilla, R., Dick, P., Donzé, F., et al. (2019). Strain and flow pathways in a shale fault zone: An in-situ test of fault seal integrity. In *Fifth International Conference on Fault and Top Seals* (pp. 8–12). Palermo, Italy: EAGE Publications.

950 Henry, W. (1803). Experiments on the quantity of gases absorbed by water, at different temperatures, and under different pressures. *Philosophical Transactions of the Royal Society of London*, **93**, 29–274.

953 Hight, D. W., Hamza, M. M., & El Sayed, A. S. (2002). Engineering characterization of the Nile Delta clays. In A. Nakase & T. Tsuchida (Eds.), *Coastal geotechnical engineering in practice* (pp. 149–162). Lisse, the Netherlands: Swets & Zeitlinger.

956 Hollingsworth, J., Jackson, J., Walker, R., & Nazari, H. (2008). Extrusion tectonics and subduction in the easter South Caspian region since 10 Ma. *Geology*, **36**(10), 763–766. https://doi.org/10.1130/G25008A.1

959 Hubbert, M. K. (1937). Theory of scale models as applied to the study of geologic structures. *Bulletin of the Geological Society of America*, **48**, 1459–1520.

961 Hustoft, S., Dugan, B., & Mienert, J. (2009). Effects of rapid sedimentation on developing the Nyegga pockmark field: Constraints from hydrological modeling and 3-D seismic data, offshore mid-Norway. *Geochemistry, Geophysics, Geosystems*, **10**(6). https://doi.org/10.1029/2009GC002409

965 Imbert, P., Geiss, B., & Fatjó de Martín, N. (2014). How to evacuate 10 km3 of mud: Saturate with gas and decrease the pressure! *Geo-Marine Letters*, **34**(2–3), 199–213. https://doi.org/10.1007/s00367-014-0357-3

968 Inan, S., Namik Yalçin, M., Guliev, I. S., Kuliev, K., & Akper Feizullayev, A. (2002). Deep petroleum occurrences in the Lower Kura Depression, South Caspian Basin, Azerbaijan: an organic geochemical and basin modeling study. *Marine and Petroleum Geology*, **14**(7–8), 731–762. https://doi.org/10.1016/s0264-8172(97)00058-5
972 Javanshir, R. J., Riley, G. W., Duppenbecker, S. J., & Abdullayev, N. (2015). Validation of lateral
973 fluid flow in an overpressured sand-shale sequence during development of Azeri-Chirag-
974 Gunashli oil field and Shah Deniz gas field: South Caspian Basin, Azerbaijan. *Marine and
975 Petroleum Geology*, *59*, 593–610. https://doi.org/10.1016/j.marpetgeo.2014.07.019

976 Jeong, S. W. (2013). Determining the viscosity and yield surface of marine sediments using
977 modified Bingham models. *Geosciences Journal*, 17(3), 241–247.
978 https://doi.org/10.1007/s12303-013-0038-7

979 Kirkham, C., Cartwright, J., Hermanrud, C., & Jebsen, C. (2017). The genesis of mud volcano
980 conduits through thick evaporite sequences. *Basin Research*, 1–20.
981 https://doi.org/10.1111/bre.12250

982 Kopf, A. J. (2002). Significance of mud volcanism. *Reviews of Geophysics*, *40*(2), 1005.
983 https://doi.org/10.1029/2000RG000093

984 Kvalstad, T. J., Andresen, L., Forsberg, C. F., Berg, K., Bryn, P., & Wangen, M. (2005). The
985 Storegga slide: Evaluation of triggering sources and slide mechanics. *Marine and Petroleum
986 Geology*, 22(1-2 SPEC. ISS.), 245–256. https://doi.org/10.1016/j.marpetgeo.2004.10.019

987 Lafuerza, S., Sultan, N., Canals, M., Lastras, G., Cattaneo, A., Frigola, J., et al. (2012). Failure
988 mechanisms of Ana Slide from geotechnical evidence, Eivissa Channel, Western
989 Mediterranean Sea. *Marine Geology*, 307–310, 1–21.
990 https://doi.org/10.1016/j.margeo.2012.02.010

991 Lawrence, G. W. M., & Cartwright, J. A. (2010). The stratigraphic and geographic distribution of
992 giant craters and remobilised sediment mounds on the mid Norway margin, and their relation
993 to long term fluid flow Slid e W. *Marine and Petroleum Geology*, 27(4), 733–747.
994 https://doi.org/10.1016/j.marpetgeo.2009.10.012

995 Locat, J., & Demers, D. (2008). Viscosity, yield stress, remolded strength, and liquidity index
996 relationships for sensitive clays. *Canadian Geotechnical Journal*, 25(4), 799–806.
997 https://doi.org/10.1139/t88-088

998 Lunne, T., Berre, T., Strandvik, S., Andersen, K. H., & Tjelta, T. I. (2001). Deepwater sample
999 disturbance due to stress relief. In *Proceedings of the 1st Annual Offshore Technology
1000 Conference* (pp. 64–85). Houston, Texas.

1001 Mainguy, M., & Longuemare, P. (2002). Coupling Fluid Flow and Rock Mechanics: Formulations
1002 of the Partial Coupling between Reservoir and Geomechanical Simulators. *Oil & Gas Science
1003 and Technology*, 57(4), 355–367.

1004 Mazzini, A., & Etiope, G. (2017). Mud volcanism: An updated review. *Earth-Science Reviews*, 168, 81–112. https://doi.org/10.1016/j.earscirev.2017.03.001

1006 Milkov, A. V. (2000). Worldwide distribution of submarine mud volcanoes and associated gas
1007 hydrates. *Marine Geology*, 167(1–2), 29–42. https://doi.org/10.1016/S0025-3227(00)00022-0

1009 Morley, C. K., von Hagke, C., Hansberry, R., Collins, A., Kanitpanyacharoen, W., & King, R.
1010 (2017). Review of major shale-dominated detachment and thrust characteristics in the
1011 diagenetic zone: Part I, meso- and macro-scopic scale. *Earth-Science Reviews*, 172(May).
1012 https://doi.org/10.1016/j.earscirev.2017.07.019
1013 Morton, A., Allen, M., Simmons, M., Spathopoulos, F., Still, J., Hinds, D., et al. (2003).
1014 Provenance patterns in a neotectonic basin: Pliocene and quaternary sediment supply to the
1015 South Caspian. Basin Research, 15(3), 321–337. https://doi.org/10.1046/j.1365-1016.2011.00208.x

1017 Mourgues, R., Bureau, D., Bodet, L., Gay, A., & Gressier, J. B. (2012). Formation of conical
1018 fractures in sedimentary basins: Experiments involving pore fluids and implications for
1019 sandstone intrusion mechanisms. Earth and Planetary Science Letters, 313–314(1), 67–78.
1020 https://doi.org/10.1016/j.epsl.2011.10.029

1021 Mukherjee, S., Talbot, C. J., & Koyi, H. A. (2010). Viscosity estimates of salt in the Hormuz and
1022 Namakdan salt diapirs, Persian Gulf. Geology Magazine, 147(4), 497–507.
1023 https://doi.org/10.1017/S001675680999077X

1024 Nermoen, A., Galland, O., Jettestuen, E., Fristad, K., Podladchikov, Y., Svensen, H., & Malthe-
1025 Srenssen, A. (2010). Experimental and analytic modeling of piercing structures. Journal
1026 of Geophysical Research: Solid Earth, 115(10), 1–15. https://doi.org/10.1029/2010JB007583

1027 Odonne, F., Imbert, P., Dupuis, M., Aliyev, A. A., Abbasov, O. R., Baloglanov, E. E., et al. (2020).
1028 Mud volcano growth by radial expansion : Examples from onshore Azerbaijan. Marine and
1029 Petroleum Geology, 112(September 2019), 104051.
1030 https://doi.org/10.1016/j.marpetgeo.2019.104051

1031 Opara, A. I. (2011). ESTIMATION OF MULTIPLE SOURCES OF OVERPRESSURES USING
1032 VERTICAL EFFECTIVE STRESS APPROACH : CASE STUDY OF THE NIGER DELTA,
1033 NIGERIA. Petroleum & Coal, 53(4), 302–314.

1034 Osborne, M. J., & Swarbrick, R. E. (1997). Mechanisms for Generating Overpressure in
1035 Sedimentary Basins : A Reevaluation 1. AAPG Bulletin, 6(6), 1023–1041.

1036 Priest, J. A., Clayton, C. R. I., & Rees, E. V. L. (2014). Potential impact of gas hydrate and its
1037 dissociation on the strength of host sediment in the Krishna-Godavari Basin. Marine and
1038 Petroleum Geology, 58(PA), 187–198. https://doi.org/10.1016/j.marpetgeo.2014.05.008

1039 Riboulot, V., Cattaneo, A., Sultan, N., Garziglia, S., Ker, S., Imbert, P., & Voisset, M. (2013). Sea-
1040 level change and free gas occurrence influencing a submarine landslide and pockmark
1041 formation and distribution in deepwater Nigeria. Earth and Planetary Science Letters, 375,
1042 78–91. https://doi.org/10.1016/j.epsl.2013.05.013

1043 Roberts, K. S., Davies, R. J., & Stewart, S. A. (2010). Structure of exhumed mud volcano feeder
1044 complexes, Azerbaijan. Basin Research, 22(4), 439–451. https://doi.org/10.1111/j.1365-1045.2011.00441.x

1046 Rudolph, M. L., Karlstrom, L., & Manga, M. (2011). A prediction of the longevity of the Lusi mud
1047 eruption, Indonesia. Earth and Planetary Science Letters, 308(1–2), 124–130.
1048 https://doi.org/10.1016/j.epsl.2011.05.037

1049 Schneider, F., Wolf, S., Faille, I., & Pot, D. (2000). A 3D Basin Model for Hydrocarbon Potential
1050 Evaluation: Application to Congo Offshore. Oil and Gas Science and Technology, 55(1), 3–13.
1051 https://doi.org/10.2516/ogst:2000001

1052 Schneider, F., Pagel, M., & Hernandez, E. (2004). Basin Modeling in a Complex Area: Example
1053 from the Eastern Venezuelan Foothills. In R. Swennen, F. Roure, & J. W. Granath (Eds.),
1054 *Deformation, fluid flow, and reservoir appraisal in foreland fold and thrust belts* (pp. 357–369). AAPG Hedberg Series. https://doi.org/10.1306/1025700H1504

1056 Sibson, R. H. (2003). Brittle-failure controls on maximum sustainable overpressure in different tectonic regimes. *American Association of Petroleum Geologists Bulletin*, 87(6), 901–908.
1057 https://doi.org/10.1306/01290300181

1059 Smith-Rouch, L. S. (2006). Oligocene-Miocene Maykop/Diatom Total Petroleum System of the South Caspian Basin Province, Azerbaijan, Iran, and Turkmenistan. *U.S. Geological Survey Bulletin*, 2201–I, 27.

1062 Stewart, S. A., & Davies, R. J. (2006). Structure and emplacement of mud volcano systems in the South Caspian Basin. *AAPG Bulletin*, 90(5), 771–786. https://doi.org/10.1306/11220505045

1064 Sultan, N., Savoye, B., Jouet, G., Leynaud, D., Cochonat, P., Henry, P., et al. (2010). Investigation of a possible submarine landslide at the Var delta front (Nice slope - SE France). *Canadian Geotechnical Journal*, 47(4), 486–496.

1067 Sultan, N., De Gennaro, V., & Puech, A. (2012). Mechanical behaviour of gas-charged marine plastic sediments. *Géotechnique*, 62(9), 751–766. https://doi.org/10.1680/geot.12.OG.002

1069 Tagiyev, M. F., Nadirov, R. S., Bagirov, E. B., & Lerche, I. (1997). Geohistory, thermal history and hydrocarbon generation history of the north-west South Caspian Basin. *Marine and Petroleum Geology*, 14(4), 363–382. https://doi.org/10.1016/S0264-8172(96)00053-0

1072 Tingay, M., Heidbach, O., Davies, R., & Swarbrick, R. (2008). Triggering of the Lusi mud eruption: Earthquake versus drilling initiation. *Geology*, 36(8), 639–642.
1074 https://doi.org/10.1130/G24697A.1

1075 Tingay, M., Manga, M., Rudolph, M. L., & Davies, R. (2017). An alternative review of facts, coincidences and past and future studies of the Lusi eruption. *Marine and Petroleum Geology*, (December). https://doi.org/10.1016/j.marpetgeo.2017.12.031

1078 Tingay, M. R. P., Rudolph, M. L., Manga, M., Davies, R. J., & Wang, C.-Y. (2015). Initiation of the Lusi mudflow disaster. *Nature Geoscience*, 8(7), 493–494.
1080 https://doi.org/10.1038/ngeo2472

1081 Torrance, J. K. (2010). Shear resistance of remoulded soils by viscometric and fall-cone methods: a comparison for the Canadian sensitive marine clays. *Canadian Geotechnical Journal*, 24(2), 318–322. https://doi.org/10.1139/t87-037

1084 Tryggvason, G. (2011). Direct Numerical Simulations of Multiphase Flows: An Introduction.

1085 Tryggvason, Gretar. (2012). *A Front-tracking/Finite-Volume Navier-Stokes Solver for Direct Numerical Simulations of Multiphase Flows*.

1087 Tryggvason, Gretar, Esmaeeli, A., Lu, J., & Biswas, S. (2006). Direct numerical simulations of gas / liquid multiphase flows. Direct numerical simulations of gas / liquid multiphase flows. *Fluid Dynamics Research*, 38, 660–681. https://doi.org/10.1016/j.fluiddyn.2005.08.006

1090 Unterseh, S., & Contet, J. (2015). Integrated geohazards assessments offshore Azerbaijan, Caspian Sea. *Offshore Technology Conference, OTC-25911* (May), 1–8.
1092 https://doi.org/10.4043/25911-MS

1093 Vincent, S. J., Davies, C. E., Richards, K., & Aliyeva, E. (2010). Contrasting Pliocene fluvial
1094 depositional systems within the rapidly subsiding South Caspian Basin; a case study of the
1095 palaeo-Volga and palaeo-Kura river systems in the Surakhany Suite, Upper Productive Series,
1096 onshore Azerbaijan. Marine and Petroleum Geology, 27(10), 2079-2106.
1097 https://doi.org/10.1016/j.marpetgeo.2010.09.007

1098 Wibberley, C. A. J., Yielding, G., & Di Toro, G. (2008). Recent advances in the understanding of fault zone internal structure: a review. *Geological Society, London, Special Publications, 299*(1), 5-33. https://doi.org/10.1144/sp299.2

1101 Wibberley, C. A. J., Gonzalez-Dunia, J., & Billon, O. (2017). Faults as barriers or channels to production-related flow: insights from case studies. *Petroleum Geoscience*, 23(1), 134-147.
1103 https://doi.org/10.1144/petgeo2016-057

1104 Woolsey, S., McCallum, M. E., & Schumm, S. A. (1975). Modeling of Diatreme Emplacement, 29-42.

1106 Wu, K., & Olson, J. E. (2015). Simultaneous Multifracture Treatments : Fully Coupled Fluid Flow and Fracture Mechanics for Horizontal Wells. In *SPE Annual Technical Conference and Exhibition* (pp. 337-346). New Orleans: Society of Petroleum Engineers.

1109 Yusifov, M., & Rabinowitz, P. D. (2004). Classification of mud volcanoes in the South Caspian Basin, offshore Azerbaijan. *Marine and Petroleum Geology*, 21, 965-975.
1111 https://doi.org/10.1016/j.marpetgeo.2004.06.002

1112 Zoporowski, A., & Miller, S. A. (2009). Modelling eruption cycles and decay of mud volcanoes. *Marine and Petroleum Geology*, 26(9), 1879-1887.
1114 https://doi.org/10.1016/j.marpetgeo.2009.03.003

1115
1116
Figure 1: Location and situation of the AMV within the South Caspian Basin background. The purple rectangle gives the location of the dense multidisciplinary dataset acquired over the Absheron anticline and integrated in the models presented in this study. A zoom of the high-resolution multi-beam bathymetry over the AMV is also displayed.

Figure 2: Updated structural model based on Green et al. (2009) work and on the fault network described in Blouin et al. (2019a). Seven layers extend along the section corresponding to different sedimentation rates, compaction laws and permeability trends (please refer to Table 1). The Layers NKG, Balakhany-Fasila, Sabunchy, Surakhany, ASF and Quaternary are named Layers 1 to 6 respectively. Numbers showed at the limits of the model correspond to limit conditions imposed for the diffusion of pore pressure and methane. Depth are given below mean sea level (MSL).

Figure 3: Results of the one-dimensional sedimentation modeling at the southern edge of the 2D geometrical model in Figure 2. On the left, vertical hydraulic conductivity versus depth trend at the end of the 5 My of sedimentation. On the right, overpressure versus depth trend at the end of the 5 My of sedimentation with corrected sedimentation rates for each layer. The top of each simulated stratigraphic unit is represented as indication using the same color code as in Figure 2 and the corresponding stratigraphic intervals are displayed in between. Depth are given below seafloor (BSF) in this figure.

Figure 4: Results of overpressure and methane diffusion modeling after 5 My of calculation considering the low permeability ASF interval. Black dotted lines are for layer limits. Layer names are reported at the right of the sections. a: overpressure (Δu) in kPa after 5 My of migration through the structural model presented on Figure 2. Overpressure migrated more rapidly through layer 4 that has a higher permeability. b: Δu/σ'v contours with values exceeding hydrofracture condition below the ASF in the north of the model, where σ'v is low. Black lines correspond to methane concentration contours. Lines are separated by 10 mM. c: Δu (kPa) vertical plot at the Absheron location (black arrow).

Figure 5: Results of overpressure and methane diffusion modeling after 5 My of calculation considering the low permeability ASF interval and faults as horizontal seals. Black dotted lines are for layer limits. Layer names are reported at the right of the sections. a: overpressure (Δu) in kPa after 5 My of migration through the structural model presented on Figure 2. Overpressure builds up along the fault network. North of the fault network overpressure is only of 18 MPa. b: Δu/σ'v contours. The highest values are now distributed south of the fault network, along the ASF, at the crest of the Absheron fold. Methane distribution is represented with black isolines, lines being separated by 10 mM. c: Δu (kPa) vertical plot at the Absheron location (black arrow).

Figure 6: Results of overpressure and methane diffusion modeling after 5 My of calculation considering the low permeability ASF interval, faults as horizontal seals and a fracture condition of 0.7. Black dotted lines are for layer limits. Layer names are reported at the right of the sections. a: overpressure (Δu) in kPa after 5 My of migration through the structural model presented on
1159 Figure 2. b: Δu/σ'v contours. Fracture occurs along the bottom edge of the ASF, south of the fault network. Methane distribution is represented with black isolines, lines being separated by 10 mM.
1160
1161 The dissolved methane distribution follows the fracture shape. c: Δu (kPa) vertical plot at the Absheron location (black arrow).
1162

1163
1164 Figure 7: Results of the simulation considering the ASF, sealing faults, a fracture condition of 0.7 and an initial methane concentration of 5550 mM after 5 My. Black dotted lines are for layer limits.
1165
1166 Layer names are reported at the right of the sections. a: overpressure (Δu) in kPa after 2 My of migration through the structural model presented on Figure 2. b: Δu/σ'v contours. Fracture occurs along the bottom edge of the ASF, south of the fault network. Methane distribution is represented with black isolines. The dissolved methane distribution follows the fracture shape and is depleted around fractures due to gas exsolution. c: degree of gas saturation (S<sub>g</sub>) calculated after fracture formation. Values as high as 1 are reached in the central part of the fracture, in an area close to the fault network. d: preconsolidation pressure (σ'<sub>p</sub>). It increases linearly with depth, but it is disturbed in the same area where gas exsolution happened reaching zero in the center of the fracture.
1167
1168
1169
1170
1171
1172
1173

1174
1175 Figure 8: Results of the simulation for μ = 10⁹ Pa.s and for the case of no overpressured mud. The
1176
1177 left column displays the evolution of density with time (left color scale) as well as the velocity
1178 vectors (m/s). The right column displays the evolution of the degree of gas saturation with time
1179
1180

1179
1180 Figure 9: Results of 1D calculations based on the case of a buoyant magma flow along a vertical dyke presented in Furbish (1997) considering a radius of 500 m corresponding to the mud source radius, compared to results obtained with 2D simulations. a: maximum velocity versus viscosity, b: minimum time for extrusion versus viscosity. Black lines with crosses correspond to the case where the initial mud overpressure is not considered, red lines are for the case with mud overpressure. Grey dots correspond to the results of the 2D simulations without mud overpressure and white dots the results of 2D simulations with mud overpressure.
1181
1182
1183
1184
1185
1186

1187
1188 Figure 10: Results of 1D calculations based on the case of a buoyant magma flow along a vertical dyke presented in Furbish (1997) considering a radius of 50 m corresponding to the conduit width, compared to results obtained with 2D simulations. a: maximum velocity versus viscosity, b: time for extrusion versus viscosity. Black lines with crosses correspond to the case where the initial mud overpressure is not considered, red lines are for the case with mud overpressure. Grey dots correspond to the results of the 2D simulations without mud overpressure and white dots the results of 2D simulations with mud overpressure.
1189
1190
1191
1192
1193
1194
1195

1195
1196 Figure 11: Formation model for the Absheron mud volcano based on in situ observations and
1197
1198 measurements, sediment analysis, laboratory testing and mud generation and remobilization
1199
1200 numerical modeling. Details of the different stages displayed in a, b, c, d, e, and f and g are in the
text. g is a zoom corresponding to black dotted rectangle in c, and shows the processes involved
in of mud remobilization towards the surface. h: legend corresponding to a, b, c, d, e and f.
1201

Table 1: Parameters used for the different diffusion simulations presented in this study. Several parameters were modified to fit observations and regional background over the different simulations. e₀ is the initial void ratio, C_c is the compression index, T the temperature, K is the hydraulic conductivity (K_h in the horizontal direction, K_v in the vertical direction), D_h is the hydraulic diffusivity, D_c the methane molecular diffusivity, Δu is the overpressure, σ'v is the vertical effective stress and e_frac the fracture void ratio.

<table><thead><tr><th></th><th></th><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 4</th></tr></thead><tbody><tr><td>nodes number*</td><td>vertically</td><td></td><td>50</td><td></td><td></td></tr><tr><td></td><td>horizontally</td><td></td><td>225</td><td></td><td></td></tr><tr><td>Seafloor temperature</td><td></td><td></td><td>5.85°C</td><td></td><td></td></tr><tr><td>grad(T)</td><td></td><td></td><td>16°C/km</td><td></td><td></td></tr><tr><td>Henry methane constant</td><td></td><td></td><td>1.5 x 10<sup>-3</sup> M.atm<sup>-1</sup></td><td></td><td></td></tr><tr><td>fluid viscosity</td><td></td><td></td><td>1.15 x 10<sup>-3</sup> Pa.s</td><td></td><td></td></tr><tr><td>methane diffusivity</td><td></td><td></td><td>1.49 x 10<sup>-7</sup> m<sup>2</sup>/s</td><td></td><td></td></tr><tr><td>e<sub>0</sub>, C<sub>c</sub>*</td><td>Layer 6</td><td></td><td>2.734; 0.159</td><td></td><td></td></tr><tr><td></td><td>Layer 5, ASF</td><td></td><td>0.500; 0.000</td><td></td><td></td></tr><tr><td></td><td>Layer 4</td><td></td><td>2.014; 0.105</td><td></td><td></td></tr><tr><td></td><td>Layer 3</td><td></td><td>1.821; 0.094</td><td></td><td></td></tr><tr><td></td><td>Layer 2</td><td></td><td>1.643; 0.080</td><td></td><td></td></tr><tr><td></td><td>Layer 1</td><td></td><td>1.659; 0.081</td><td></td><td></td></tr><tr><td>a*, b*, K<sub>h</sub>/K<sub>v</sub>*</td><td>Layer 6</td><td></td><td>3.064; -22.866; 2200.0</td><td></td><td></td></tr><tr><td>K<sub>h</sub>=exp(ae+b)</td><td>Layer 5, ASF</td><td></td><td>0.000; -25.328; 15.0</td><td></td><td></td></tr><tr><td>(m/s)</td><td>Layer 4</td><td></td><td>3.064; -17.666; 15.0</td><td></td><td></td></tr><tr><td></td><td>Layer 3</td><td></td><td>3.064; -22.150; 5.0</td><td></td><td></td></tr><tr><td></td><td>Layer 2</td><td></td><td>3.064; -21.803; 2.0</td><td></td><td></td></tr><tr><td></td><td>Layer 1</td><td></td><td>3.064; -22.333; 1.0</td><td></td><td></td></tr><tr><td>e<sub>0</sub>, C<sub>c</sub>*</td><td>faults</td><td></td><td>3.5; 0.0</td><td></td><td></td></tr><tr><td>D<sub>h</sub> (m<sup>2</sup>/s)</td><td></td><td>5.5 x 10<sup>-6</sup></td><td>5.5 x 10<sup>-8</sup></td><td>5.5 x 10<sup>-8</sup></td><td>5.5 x 10<sup>-8</sup></td></tr><tr><td>D<sub>c</sub> (m<sup>2</sup>/s)</td><td></td><td></td><td>5.5 x 10<sup>-6</sup> m<sup>2</sup>/s</td><td></td><td></td></tr><tr><td>fracture condition: Δu/σ'<sub>v</sub>*</td><td></td><td>0.9</td><td>0.9</td><td>0.7</td><td>0.7</td></tr><tr><td>e<sub>frac</sub>*</td><td></td><td></td><td>1.0</td><td></td><td></td></tr><tr><td>Start of overpressure diffusion (Ma)</td><td></td><td>1 My</td><td>1 My</td><td>1 My</td><td>3 My</td></tr><tr><td>Initial methane concentration (mM)</td><td></td><td>55.5</td><td>55.5</td><td>55.5</td><td>5550</td></tr><tr><td colspan="6">*: Dimensionless values</td></tr></tbody></table>

1208

Table 2: Parameters used in the different simulations completed during the study. The varying parameters are the sole sediment viscosity and the fracture length.

<table><thead><tr><th rowspan="2">Mud condition</th><th rowspan="2">viscosity (Pa.s)</th><th rowspan="2">conduit length (m)</th><th colspan="3">dimensions</th><th colspan="2">mud source</th><th colspan="4">Mass-densities (kg/m<sup>3</sup>)</th></tr><tr><th>length</th><th>depth</th><th>dimension</th><th>position</th><th>overpressure</th><th>bottom</th><th>top</th><th>mud</th><th>Solid grains</th></tr></thead><tbody><tr><td rowspan="3">No overpressured mud</td><td rowspan="3">10<sup>8</sup></td><td>0</td><td>10 km</td><td>3900</td><td>radius 500</td><td>x =</td><td>0 MPa</td><td rowspan="3">2100</td><td rowspan="3">1900</td><td rowspan="3">1900</td><td rowspan="3">2650</td></tr><tr><td>1000</td><td>100 nodes</td><td>m</td><td>m</td><td>5000 m</td><td></td></tr><tr><td>2000</td><td>100 nodes</td><td>m</td><td>m</td><td>3900 m</td><td></td></tr><tr><td></td><td>10<sup>9</sup></td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>1000</td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><tbody><tr><td></td><td>10<sup>10</sup></td><td>2000</td></tr><tr><td></td><td></td><td>0</td></tr><tr><td></td><td></td><td>1000</td></tr><tr><td></td><td></td><td>2000</td></tr><tr><td rowspan="3">Overpressured mud</td><td>10<sup>8</sup></td><td>0</td></tr><tr><td></td><td>1000</td></tr><tr><td></td><td>2000</td></tr><tr><td rowspan="3">10<sup>9</sup></td><td></td><td>0</td></tr><tr><td></td><td>1000</td></tr><tr><td></td><td>2000</td></tr><tr><td rowspan="3">10<sup>10</sup></td><td></td><td>0</td></tr><tr><td></td><td>1000</td></tr><tr><td></td><td>2000</td></tr></tbody></table>

35 MPa

1211

<table><caption>Table 3: Synthesis of the main results obtained from the different simulations computed in this study. The initiation time of gassy-mud ascent, the time needed to reach the seafloor, the maximum velocity and the final diameter of the gassy-mud column are displayed for the two initial conditions concerning mud pressures.</caption><thead><tr><th>Mud condition</th><th>μ (Pa.s)</th><th>Conduit length (m)</th><th>Initiation (s)</th><th>Extrusion (s)</th><th>Max. velocity (m/s)</th><th>Final mud column diameter (m)</th></tr></thead><tbody><tr><td rowspan="6">no overpressured mud</td><td rowspan="3">10<sup>8</sup></td><td>0</td><td>350</td><td>600</td><td>3.2</td><td>1000</td></tr><tr><td>1000</td><td>300</td><td>750</td><td>1.2</td><td>250</td></tr><tr><td>2000</td><td>300</td><td>750</td><td>0.5</td><td>100</td></tr><tr><td rowspan="3">10<sup>9</sup></td><td>0</td><td>3500</td><td>6000</td><td>0.25</td><td>500</td></tr><tr><td>1000</td><td>3000</td><td>7000</td><td>0.12</td><td>250</td></tr><tr><td>2000</td><td>3000</td><td>7500</td><td>0.05</td><td>100</td></tr><tr><td rowspan="6">overpressured mud</td><td rowspan="3">10<sup>10</sup></td><td>0</td><td colspan="4">Long computation time – no output data</td></tr><tr><td>1000</td><td colspan="4"></td></tr><tr><td>2000</td><td colspan="4"></td></tr><tr><td rowspan="3">10<sup>8</sup></td><td>0</td><td>200</td><td>400</td><td>3.5</td><td>750</td></tr><tr><td>1000</td><td>200</td><td>500</td><td>4.2</td><td>400</td></tr><tr><td>2000</td><td>200</td><td>500</td><td>3.2</td><td>100</td></tr><tr><td rowspan="3">10<sup>9</sup></td><td>0</td><td>1750</td><td>4000</td><td>0.35</td><td>750</td></tr><tr><td>1000</td><td>1500</td><td>4500</td><td>0.36</td><td>400</td></tr><tr><td>2000</td><td>2000</td><td>4750</td><td>0.32</td><td>100</td></tr><tr><td rowspan="3">10<sup>10</sup></td><td>0</td><td colspan="4">Long computation time – no output data</td></tr><tr><td>1000</td><td colspan="4"></td></tr><tr><td>2000</td><td colspan="4"></td></tr></tbody></table>

1216
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>
Absheron

a S

Depth (m)

b

Depth (m)

c

Depth (m)

d

Depth (m)

N

L6
L5
L4
L3
L2
L1

Maykop

L6
L5
L4
L3
L2
L1

Maykop

L6
L5
L4
L3
L2
L1

Maykop

L6
L5
L4
L3
L2
L1

Maykop

<figure><img src="image_8.png" /></figure>

500

<figure><img src="image_9.png" /></figure>

Methane
concentration

Horizontal distance (m)

Δu (kPa)

<figure><img src="image_10.png" /></figure>

Δu/σ'<sub>v</sub>

<figure><img src="image_11.png" /></figure>

S<sub>g</sub>

<figure><img src="image_12.png" /></figure>

σ'<sub>p</sub> (kPa)

<figure><img src="image_13.png" /></figure>

Accepted Article
<figure><img src="image_14.png" /></figure>

0 sec

0 sec

<figure><img src="image_15.png" /></figure>

3750 sec

3750 sec

<figure><img src="image_16.png" /></figure>

4500 sec

4500 sec

<figure><img src="image_17.png" /></figure>

6000 sec

6000 sec

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
a

<figure><img src="image_20.png" /></figure>

b

<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /></figure>

a

Time for extrusion (years)

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /></figure>

g

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

h