Methane bubble growth in fine-grained muddy aquatic sediment: Insight from modeling

Earth and Planetary Science Letters 377–378 (2013) 336–346

Contents lists available at SciVerse ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Methane bubble growth in fine-grained muddy aquatic sediment: Insight from modeling Regina Katsman a,∗ , Ilia Ostrovsky b , Yizhaq Makovsky a a b

The Dr. Moses Strauss Department of Marine Geosciences, Faculty of Science and Science Education, The University of Haifa, Haifa, Mount Carmel 31905, Israel Yigal Alon Kinneret Limnological Laboratory, Israel Oceanographic and Limnological Research, P.O.B. 447, Migdal 14950, Israel

a r t i c l e

i n f o

Article history: Received 9 November 2012 Received in revised form 2 May 2013 Accepted 7 July 2013 Available online 12 August 2013 Editor: L. Stixrude Keywords: bubble methane muddy sediment bubble growth model fracture mechanics solute transport

a b s t r a c t Methane (CH4 ) is the most abundant hydrocarbon and one of the most important greenhouse gases in the atmosphere. CH4 bubble growth and migration within muddy aquatic sediments are closely associated with sediment fracturing. In this paper we present the modeling of buoyancy-driven CH4 bubble growth in fine-grained muddy aquatic sediment prior to the beginning of its rise. We designed a coupled mechanical/reaction-transport numerical model that enables a differential fracturing over the bubble front (as it occurs in nature), when the fracturing increment stays constant at the bubble head and subsides towards bubble tail during bubble growth. We show that this differential fracturing over the bubble front controls the bubble shape and size temporal evolution, and is significantly affected by the critical stress intensity factor of the muddy sediment. The intercalated stages of elastic expansion and fracturing during the bubble growth shorten with time as the bubble approaches its terminal size (prior to its ascent). Our simulations reveal a high asymmetry in the bubble shape growing with time, with respect to its initial symmetric penny-shaped configuration. It was found that the bubble grows allometrically, while the importance of the bubble surface area growth with time. We also confirmed the earlier predictions about the ”inverted tear-drop” bubble cross-section just prior to the beginning of its rise. Modeling of the terminal bubble characteristics will permit prediction of the delivery of gaseous methane from the sediment to the atmosphere via the water column. © 2013 Elsevier B.V. All rights reserved.

1. Introduction

Methane (CH4 ) is the most abundant hydrocarbon and one of the most important greenhouse gases in the atmosphere. Over the last century the CH4 concentration has risen by 1% per year (Rowland, 1985). Despite their importance, CH4 fluxes from the aquatic sediments recently reported in the literature ranged over an order of magnitude (Soumis et al., 2005; St. Louis et al., 2000) indicating that they have not been properly quantified (Del Sontro et al., 2010). CH4 emission from aquatic systems to the atmosphere is usually dominated by gas ebullition (Ostrovsky et al., 2008; Del Sontro et al., 2010): In shallow lakes up to 98% of CH4 release originated from bubbles while only 2% came from dissolved CH4 (Keller and Stallard, 1994; Del Sontro et al., 2010). Formation of methane bubbles and their transport within the sediments is a subject of recent investigations.

*

Corresponding author. Tel.: +972 4 8288979; fax: +972 4 8288267. E-mail addresses: [email protected] (R. Katsman), [email protected] (I. Ostrovsky), [email protected] (Y. Makovsky). 0012-821X/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.epsl.2013.07.011

1.1. The growth and migration of bubbles within muddy sediments The perceived pliability of soft muddy sediments and the observed fluidization patterns (e.g. gravity flow), erroneously suggest that such sediments could act fluidly or plastically in response to stress induced by growing bubble (Wheeler, 1988). However, recent laboratory simulations have shown that these sediments respond mechanically as fracturing elastic solid (Best et al., 2004; Boudreau et al., 2005; Barry et al., 2010). The importance of grain size in determining the behavior of gassy sediments has recently been demonstrated by Jain and Juanes (2009), and Choi et al. (2011). They suggested that gas migration in fine-grained (muddy) sediments is governed by a fracture-dominated regime due to the large capillary-entry pressure precluding gas from entering pore throats without breaking the inter-granular bonds, while in coarsegrained (sandy) sediments it occurs by capillary invasion through the sediment framework. Exploring the bubble shape, Anderson et al. (1998) observed that bubbles in shallow muddy sediment were often non-spherical, with eccentricity increasing with their volume. Imaging bubble in gassy sediments, Best et al. (2004) demonstrated that in contrast to the small (2 mm in diameter) sub-spherical bubbles in silty sands, most bubbles observed in the clayey silts appeared as low aspect

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

337

ratio cavities, up to ∼40 mm long, with their longest axes aligned in the sub-vertical plane. Van Kessel and van Kesteren (2002), Winterwerp and van Kesteren (2004), Johnson et al. (2002) demonstrated experimentally that a bubble initially grows quickly to entirely fill the pore space it occupies, and then starts deforming the surrounding sediment skeleton. After numerous observations of crack initiation and propagation in soft sediments, it was suggested that these cracks can be described by linear elastic fracture mechanics (LEFM), the theory applicable to cracks in which the region of plastic strain at their tips is small relative to the crack size (Lawn and Wilshaw, 1975; Broek, 1986). 1.2. Recent modeling approaches Recent modeling efforts addressed the issue of methane dynamics in sediments through several alternative perspectives. 1.2.1. Diagenetic reaction-transport models Despite the importance of exploring gaseous methane dynamics, models incorporating an explicit gas phase representation usually use a multiphase fluid dynamics approach that is better justified for coarse-grained soils (Oldenburg et al., 2010) or fractured aquifers (Rubin et al., 2008). Martens et al. (1998) incorporated a gas phase into a diagenetic model, balancing in a steady-state sedimentation and methane production from organic matter against ebullition of gaseous methane, CH4 ( g ). Haeckel et al. (2004) represented CH4 ( g ) as a source term for dissolved methane, implying that CH4 ( g ) is not transported explicitly through the sediment, while Davie and Buffett (2001) assumed that gas transport follows burial. Alternatively, Haeckel et al. (2007) considered gas phase transport explicitly through tube structures. Further, Dale et al. (2008) added a mass-conservation equation for CH4 ( g ), assuming that exchange of methane between the dissolved and gaseous phases is proportional to the departure from the local solubility (Duan et al., 1992). Mogollon et al. (2009, 2011) derived a separate 1D mass and momentum conservation equations for the solid, aqueous, and gas phases coexisting within a common control volume, in addition to conservation of the individual species. 1.2.2. Mechanical models Microscopic Discrete-Element Model coupling two-phase (gas– brine) flow with sediment mechanics was implemented in Jain and Juanes (2009) allowing sediment fracture by an advancing gas phase, caused by localized breaking of cohesion between adjacent sediment grains. A steady state (Gardiner et al., 2003a) and transient (Algar and Boudreau, 2009, 2010) reaction–diffusion model describing supply of dissolved methane, CH4 (aq), to the growing bubble was combined with principles of LEFM to simulate sediment elastic expansion, followed by uniform fracturing at the bubble front during its growth. A key assumption involved in the model of Gardiner et al. (2003a) is that bubbles grow slow enough to enable the solute concentration field next to the bubble be readily adjusted. This assumption allowed Gardiner et al. (2003b) to find an analytical solution to the diffusion equation in oblate-spheroidal coordinates. However, fracture event results in sudden increase in bubble volume and drop in internal gas pressure (Johnson et al., 2002), implying that transient reaction–diffusion equation must be solved in combination with LEFM, as implemented by Algar and Boudreau (2009, 2010). When buoyancy was added to the model, an initial rise (propagation) of a large (mature) bubble was simulated by Algar et al. (2011b), neglecting mass transfer between the bubble and sediment. Fracturing in the model was permitted to occur at the bubble head only.

Fig. 1. Sediment cell is presented as a block cut by a symmetry plane (only half a block is modeled). The initial bubble seed is presented as a small penny-shaped crack placed on the symmetry plane. Only one of the two bubble free surfaces is modeled. Points B and A are the head and tail points on the bubble front respectively. Point P prescribed on the bubble surface is used for the calculations of Stress Intensity Factor at the adjacent point F at the bubble front (see Section 2.1.3).

Coupling mechanical and reaction-transport processes is the issue addressed by a structural diagenesis (Laubach et al., 2010). This coupling is greatly important for general understanding of bubble growth and migration in fine-grained muddy sediments. In general, such a feedback can help to obtain a scientific knowledge about the low-temperature realm of sedimentary basins that is of great intrinsic and practical interest. However, despite the undoubted progress in development of such models, the accurate reproduction of the processes governing bubble dynamics is still scarce. In this paper we present a model describing the process of the buoyancy-driven bubble growth in muddy sediments prior to beginning of its ascent toward the sediment-water interface. For the first time we simulate the differential fracturing over the bubble front, as it occurs in nature, and show that the fracturing controls the bubble shape and size evolution. The model has a large potential for simulating buoyancy-driven upward migration of bubble, and predicting the delivery of gaseous methane from sediments to the atmosphere. 2. Methods 2.1. The model To model bubble growth within soft sediment we used approach of bubble growth within an elementary cell that was applied by Proussevitch et al. (1993), Proussevitch and Sahagian (1996), Navon et al. (1998), Favelukis (2004) for fluids, and by Algar and Boudreau (2009, 2010) for muddy sediments. Our model extends the above approaches by including a fully-coupled numerical treatment of transport and mechanical components. It especially concentrates on the precise modeling of the differential fracturing of muddy sediments at the growing bubble front, which is carried out accordingly to the principles of fracture mechanics. The modeling setup is depicted in Fig. 1. The sediment cell is presented as a block cut by a symmetry plane (only half of the block is modeled). The initial bubble seed is presented as a small pennyshaped crack with one of its free surfaces located on the symmetry plane (Algar et al., 2011b). 2.1.1. Modeling of the region outside the bubble Transport of dissolved methane and solid mechanics are modeled in the bulk of the sediment cell outside the bubble (Fig. 1).

338

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

Solute transport The Solute Conservation Equation in the sediment is modeled for the dissolved methane, CH4 (aq), to capture the dynamics of the solute supply to the growing bubble (e.g. Lasaga, 1997). Only diffusive transport within the elementary cell is considered at this stage:

  ∂ C CH4 (aq) + ∇ · − D ∇ C CH4 (aq) = S CH4 (aq) ∂t

(1)

where C CH4 (aq) is aqueous methane concentration, D is tortuositycorrected methane diffusion coefficient (Boudreau, 1997), S CH4 (aq) is dissolved methane production/consumption rate (a source strength), and t is time. Boundary conditions: Owing to the relatively small dimensions of the computational domain located far below SMTZ, the difference in a dissolved CH4 (aq) concentrations between the boundaries at the prescribed location is small (see Martens et al., 1998) and may be neglected (as in Algar and Boudreau, 2009). As a result, a uniform solute concentration is assumed at all the boundaries of the domain, C CH4 (aq) = C 0 , as in Proussevitch et al. (1993), Proussevitch and Sahagian (1996), Navon et al. (1998), except at the symmetry plane. No flux boundary condition is specified at  is n(− D ∇ C CH4 (aq)) = 0 (Fig. 1), where n the symmetry plane, − outward normal to the surface, D is methane diffusion coefficient. Solid mechanics Algar et al. (2011b) suggested that the behavior of the muddy sediment can be described by a Kelvin–Voigt model. At the longterm timescale it is fully described by LEFM, and at the short timescale it is captured by the added viscosity. This choice was motivated by the overall mechanical stability of bubbles within muddy sediment lacking long-term deformation. If the time required for the information about stress changes to travel to the crack tail (τ f = H /U R , where H is bubble height and U R is Raleigh velocity of the material, Algar et al., 2011b) is larger than the retardation time of Kelvin–Voigt model (τ v = μ/ E) (Malvern, 1969), then the model becomes purely elastic. It is demonstrated in Section 2.3 that implementation of the LEFM approach (neglecting the viscous component) is justified in the current model. Moreover, as the main focus of this study is to model evolution of size and shape of the growing bubble, the time-dependent aspect of the viscous material response reveals no factual importance. Thus, the developed solid mechanics part of the model includes two major components: elasticity and fracturing. (a) Linear elasticity is modeled with the Force Equilibrium Equation, simulating elastic deformations caused by the growing pressure within a bubble and by remote loading:

 ∂2 w  ) = F g ρ 2 − ∇ · (c ∇ w ∂t

(2)

 is a vector of elastic displacewhere c is a stiffness tensor, w ment, F g is gravity load, and ρ is material bulk density. Boundary conditions: Vertical compressive stress, F , is applied onto the upper horizontal face of the computational domain, σ · n = − F , where σ is a stress tensor. Zero normal displacement is prescribed onto all vertical faces of the domain, · w  = 0, including a symmetry plane. The bottom face is fixed, n  = 0. w (b) Fracturing is described in Section 2.1.3 below. 2.1.2. Modeling of bubble Bubble nucleation (the appearance of the initial crack) is ignored in the model following the assumptions of the LEFM (Lawn and Wilshaw, 1975).

Conservation of methane for the growing bubble is calculated by integrating the dissolved methane flux over the bubble surface, α (e.g. Algar and Boudreau, 2009):

∂(C CH4 ( g ) V b ) = ∂t







 · D φ∇ C CH4 (aq) dα n

(3)

α

where C CH4 ( g ) and C CH4 (aq) are gaseous and aqueous methane  is outward concentrations, respectively, V b is bubble volume, n normal at the bubble surface α , φ is effective porosity, D is tortuosity-corrected methane diffusion coefficient. Pore-water dissolved methane concentration, C CH4 (aq), at the bubble surface α is in equilibrium with gaseous methane concentration inside the bubble, C CH4 ( g ), according to Henry Law:

C CH4 ( g ) = C CH4 (aq)(α , t ) · k H

(4)

where k H is dimensionless Henry Law constant (Gardiner et al., 2003a, 2003b; Algar and Boudreau, 2010). At any stage of the elastic loading the bubble volume is calculated by integrating elastic  in the normal direction, n over the bubble surdisplacements, w, face:



·w  dα n

Vb =

(5)

α

Methane is treated as an ideal gas, and its concentration, C CH4 ( g ), derived from Eq. (4) along with capillary pressure term is used to calculate the pressure within the bubble:

P = C CH4 ( g ) R g T + 2γ /rm

(6)

 is outward normal at where R g is gas constant, T is temperature, n the bubble surface, γ is surface tension. rm is the geometric mean radius of bubble surface curvature, r2 = ( r1 + r1 ) (Pinder and Gray, m 1 2 2008). r1 and r2 are the radii of curvature of two orthogonal curves prescribed on the surface (Fig. 1, dashed lines) and calculated at each moment of bubble growth. This capillary pressure term is associated with a thin coating of water along the bubble surfaces between the gas and the sediment, if present, rather than with water at pore throats in fine-grained muddy sediments (Boudreau, 2012). Gaseous methane pressure, P , calculated in Eq. (6) is applied as the boundary condition on the bubble/sediment interface (i.e. on the bubble free surface, Fig. 1):

σ · n = − P

(7)

2.1.3. Modeling of fracturing Cracking induced by pressure growth within the bubble is modeled using principles of Linear Elastic Fracture Mechanics. LEFM postulates that the behavior of crack is determined solely by the value of the stress intensity factor (SIF), K , at its front (Lawn and Wilshaw, 1975). In contrast to cracks in 2D, even for the simple elliptical 3D crack SIF, varies along the crack front even under onedimensional symmetric remote extension (Kassir and Sih, 1975). Irwin (1957) derived a relationship between the strain energy release rate and the SIF, G ∝ K 2 , providing a link between the crack-tip stress field and the energy-balance criterion for crack growth. SIF at the crack (bubble) front is evaluated in this paper using Displacement/Stress Extrapolation Method, from the known displacement in the vicinity of the crack front, utilizing the general analytical solution for the crack problem (Broek, 1986; Kassir and Sih, 1975). Here, crack opening displacements are evaluated on the crack’s surface (Fig. 1) using a one point methodology (e.g. Couroneau and Royer, 2000; Guo and Shen, 2003). In this

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

339

method point P is prescribed on the crack surface in the vicinity of point F at the crack front with distance r from it (Fig. 1). Mode I SIF at each point F of the front is calculated using:

KI =

E 4(1 − ν 2 )



π 2r

2w nP

(8)

 P at point P , Here w nP is the projection of the displacement w  normal to the crack’s surface, α , E is Young’s in the direction n modulus, and ν is Poisson ratio. It can easily be verified that under the loading induced by the growing bubble K II and K III are neglected. The crack propagation direction is found using the classical Strain Energy Density Criterion (Sih, 1974). In the general case of loading strain energy density factor, is given by: Se(θ) = a11 K I2 + 2 2a12 K I K II + a22 K II2 + a33 K III , where the coefficients ai j are defined by Sih (1974) as a function of θ . Here θ is an out of plane angle − → measured from the direction normal to the crack front, nn, in the local crack coordinate system (Sih, 1974, Fig. 1; Citarella and Cricri, 2010, Fig. 5), i.e. θ = ±π coincide with crack surfaces, α . The theory is based on three hypotheses: (1) The direction of crack growth at any point on the cracks front is towards the region with minimum strain energy density factor, Se(θ) = Semin (θ), as compared with other directions on the same spherical surface surrounding the point; (2) Crack extension occurs when Se(θ) in the region reaches a critical value, Sec (θ); (3) The variable length, a, of the initial crack increment is assumed to be proportional to Semin (θ), such that Semin (θ)/ a = (dW /dV )cr remains constant along the new crack’s front, where (dW /dV )cr is a critical strain energy per unit volume, a material property (see Fig. 6 in Kassir and Sih, 1975). In our case with K II = K III = 0, the strain energy density factor is Se(θ) = a11 K I2 , with coefficient a11 = (1 + cos(θ))/16η · [2 · (1 − 2ν ) + 1 − cos(θ)] (Sih, 1974). Here η is the shear modulus, ν is Poisson ratio. By minimizing the latter equation with respect to θ one can obtain that min(a11 ) is attained at θ = 0, which ensures Mode I crack propagation. According to the point (3) above, the distribution of crack increments along the crack front is calculated as

a = amax ( K I / K I max )2

(9)

(Citarella and Cricri, 2010), where K I max is maximal SIF at the crack front that exceeds a prescribed critical SIF K I max  K Ic to ensure the fracturing. The maximum crack increment amax is prescribed here to be much smaller than the size of a plastic zone at the crack tip used in Algar and Boudreau (2009). Following the methodology of fracture mechanics simulations, amax is prescribed as the size of the mesh element that is adjacent to the crack front (Citarella and Cricri, 2010; Tabiei and Wu, 2003). Once the criterion for the crack growth is met, the crack propagates in a prescribed direction by a given local increment a calculated at its front. This process is modeled numerically using mesh deformation. The computational mesh is moved by a prescribed variable displacement a at the bubble front (Fig. 1), while a new mesh is not immediately generated at each time step. Instead, the model perturbs the mesh nodes within the computational domain so that they conform to the deformed boundary (bubble free surface). The mesh in the domain is smoothed using the Laplace smoothing method. Only when mesh quality becomes poor is remeshing of the geometry conducted. All the equations specified outside and inside the bubble are simultaneously solved in a fully coupled manner. The model is built using the Comsol Multiphysics simulation environment.

Fig. 2. Model verification. (a) Cross-section of penny-shaped crack cut at X = 0 (Fig. 1). R is crack semi-axis. Left panel: Initial crack configuration when the internal pressure, P , equals to the remote external stress, σ y ( P = σ y ). Right panel: Inflated crack under P > σ y . Crack opening displacement (COD) is measured at distance z from the crack center. (b) Distribution of COD/2 calculated by numerical and analytical methods and depicted as a function of coordinate z scaled with crack semi-axis, R. A good agreement between the solutions is found.

2.2. Model verification The elastic part of the model was verified against the analytical solution for the Crack Opening Displacement (COD) for a pennyshaped (3D) crack (Sneddon, 1946; Kassir and Sih, 1975):

√

COD =

8(1 − ν 2 )( P − σ y )

πE

R 2 − z2

(10)

where R is crack semi-axis, z is the vertical spatial coordinate with origin at the crack center, P is an internal bubble pressure, and σ y is the remote stress component acting in the direction perpendicular to crack surfaces (Fig. 2(a)), ν is Poisson ratio, and E is Young’s modulus. Under conditions of uniaxial strain (e.g. Turcotte and Shubert, 1982) prescribed in this study σ y is calculated as σ y = 1−ν ν σz , where σz is vertical remote stress component. Verification of the model was conducted on an initial model configuration (Fig. 1) with geometrical and material parameters presented in Section 2.3, with vertical remote stress σz = 1.5 kPa, and bubble internal pressure P = 3 kPa. The distribution of COD/2 with coordinate z, scaled with crack semi-axis R, is shown in Fig. 2(b) for both numerical and the analytical solutions. A good agreement between the solutions is demonstrated.

340

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

2.3. Simulation conditions A computational domain is prescribed with equal height and length, h = 25 cm, and a smaller width b = 8.3 cm (Fig. 1). A computation starts with small flat penny-shaped bubble with initial radius R = 4 mm to ensure a minimal difference in SIF between the tail (A) and the head (B) points at the crack front (Fig. 1), when inflated. Bubble growth before this stage (when R < 4 mm) is ignored as beyond the pore scale the small bubble retains a penny-shaped configuration (Boudreau, 2012). Ambient conditions for the calculations are taken as those of the NRL site in Eckernförde Bay (26 m water depth) at 1 m sediment depth for the summer conditions:

• Dissolved CH4 (aq) (saturated) methane concentration is C 0 = 6.5 mM = 0.1 kg m−3 (Martens et al., 1998). This value is applied as boundary and initial conditions for the calculations.

• Source strength is S CH4 (aq) = 2.56 e−10 kg m−3 s−1 (Martens et al., 1998).

• Diffusion coefficient (tortuosity-corrected) is approximated as D = 5 · 10−10 m2 s−1 (Algar and Boudreau, 2009). • Sediment bulk density is taken as ρ = 1240 kg m−3 (Silva and Brandes, 1998; Anderson et al., 1998; Richardson and Briggs, 1996), and effective overburden stress as σz = 1.5 kPa (Silva and Brandes, 1998). • Effective porosity is taken as φ ∼ = 0.5 based on the following considerations. Depending on the mineralogy of the finegrained clays, the immobile double layer of water adsorbed on the clay surface and dead-end pores usually eliminate a portion of the interconnected pore space where solute movement can occur. Thus, effective porosity of clays may range between 25% and 100% of the total porosity (Seevee, 2010), while in some clays it may approach zero. Based on the available data (Martens et al., 1998; Mogollon et al., 2009) total porosity at the given location is estimated as φ T ∼ = 0.8, and effective porosity is assumed to be 60% of this value, i.e. φ ∼ = 0.5. • Young’s modulus is taken as E = 1 · 105 Pa (Algar and Boudreau, 2010), and critical stress intensity factor is K c = 200 Pa · m1/2 . • Poisson ratio is ν = 0.45 (Dorgan et al., 2007; L’Esperance et al., 2012). The time scale of the elastic response is defined as τ f = H /U R (Algar et al., 2011b), where H is bubble height and U R is Raleigh wave velocity of the material (Gross and Seelig, 2006). In our case maximal bubble height is H = 13 cm (see Section 5 and Fig. 7), and U R can be approximated as U R = 0.955U E for the current type of the sediment (Algar et al., 2011b), where U E is shear wave velocity for muddy sediments approximated for the prescribed location in the Eckernförde Bay as U E ∼ = 12 m/s (Richardson and Briggs, 1996; Davis et al., 1996), and confirmed by theoretical calculations by Buckingham (2005). These values result in τ f = 10−2 s. The retardation time in the Kelvin–Voigt model is τ v = μ/ E (Malvern, 1969), where μ is dynamic viscosity and E is Young’s modulus. Considering kinematic viscosity that varies between 20 and 2000 cm2 /s for marine mud (Jiang and LeBlond, 1993; Mei and Liu, 1987, and references therein), which contrasts to viscosities of mudslides or debris flows, and E = 1.e5 Pa, the retardation time is always smaller than the time scale of the elastic response (i.e. τ v < τ f ). This justifies the validity of the LEFM model within the considered range of the viscosities. 3. Results In this section we analyze the evolution of bubble shape and size under buoyant loading in muddy sediment prior to its rise.

Fig. 3. Distribution of SIF at the bubble front (to the right from the vertical symmetry line connecting points A and B in Fig. 1). For the small initial penny-shaped bubble with radius R = 0.004 m the difference in SIF between points A and B reaches K IAB = 2 Pa · m1/2 (the red line). It is depicted at the moment when SIF at B approaches a critical value of K Ic = 200 Pa · m1/2 prior to the fracturing initiation (t = 444 s). K IAB grows with time as fracturing proceeds (the blue line). At t = 1349 s it approaches K IAB = 8 Pa · m1/2 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 demonstrates distribution of SIF at the bubble front (to the right from the vertical symmetry line connecting points A (bubble tail) and B (bubble head), Fig. 1). For the small initial penny-shaped bubble with radius R = 0.004 m the difference in SIF between points A and B is minimal and reaches K IAB = 2 Pa m1/2 (the red line) at the moment when SIF approaches a critical value of K Ic = 200 Pa m1/2 at the crack head (B) prior to fracturing initiation (t = 444 s). K IAB grows with time as fracturing proceeds (the blue line). For instance, at t = 1349 s the difference approaches K IAB = 8 Pa m1/2 . Fig. 4 demonstrates the evolution of pressure in the bubble with time. During the first 444 s, initial circular bubble grows exclusively by an elastic expansion. When the critical value of the SIF, K Ic = 200 Pa m1/2 , is attained at the bubble head, fracturing starts at the entire bubble front. Each fracturing event is accompanied by an instantaneous drop in pressure within the bubble, followed by a short-term gradual rise in pressure induced by solute transport to the growing bubble. A general decrease in the internal bubble pressure with time is also observed. Zooming into the details of the internal bubble pressure evolution (Fig. 5) reveals that both pressure drop increment and the duration of the elastic expansion increments are different at the beginning and at the end of the bubble growth. Specifically, the pressure drop increment in the beginning of the bubble growth comprises about 20 Pa, while about 1 Pa at the end. This difference can be attributed to the fact that SIF at the bubble head approaches a critical value as the bubble grows (as it follows from the principles of fracture mechanics, e.g. Lawn and Wilshaw, 1975). In addition, solute supply step extends for about 2.3 s in the beginning of the simulations and about 0.025 s in the end, thus indicating a more effective solute supply via the bubble surface for the larger bubble. The capillary pressure (Fig. 6) is found to be negligible compared to the pressure developed within the bubble (Fig. 4), even prior to the initial fracturing event (t = 444 s). It subsides even more with time as fracturing proceeds and the bubble grows. The evolution of crack increment, a, at the bubble front is shown in Fig. 7 as a function of the vertical coordinate z. This

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

Fig. 4. The evolution of pressure in the bubble with time. During the first 444 s, initial circular bubble grows exclusively by an elastic expansion. When the critical value of the SIF, K Ic = 200 Pa m1/2 , is attained at the bubble head, fracturing starts on the entire bubble front. Each fracturing event is accompanied by an instantaneous drop in pressure within the bubble, followed by a short-term gradual rise in pressure induced by solute transport to the growing bubble. A general decrease in the internal bubble pressure with time is also observed.

Fig. 5. Zoom into the details of internal pressure evolution at the beginning and at the end of the bubble growth. The pressure drop increment in the beginning of the bubble growth comprises about 20 Pa, while about 1 Pa at the end. Solute supply step extends for about 2.3 s in the beginning of the simulations and about 0.025 s in the end. See the text for the discussion on the matter.

parameter controls bubble shape and size evolution. The time period from the first fracturing event (t = 444 s) until the moment when a at the tail of the bubble drops to zero (t = 3524 s) is presented. The top curve corresponds to the initial (t = 444 s) distribution of the crack increment at the front when z A = −0.004 m z B = 0.004 m (vertical coordinates of the bubble tail and head). The bottom curve corresponds to final moment of the fracturing, when a at the bubble tail drops to zero (t = 3524 s) when z A = −0.039 m, z B = 0.09 m. Fig. 7 demonstrates that crack propagation increment always remains constant at the bubble head and subsides at the bubble tail. After t = 3524 s bubble never grows at its tail. Fig. 8 demonstrates the initial (t = 444 s) and final (t = 3524 s) bubble shape and size in the front view (panel a, X–Z plane in Fig. 1), and in the cross-section (panel b, Z –Y plane in Fig. 1). Both panels show that bubble grew significantly and became asymmetrical with respect to the initial penny-shaped configuration. The bubble opening (panel b) is depicted by the normal displacement component, w n , perpendicular to the bubble surface. This cross-

341

Fig. 6. The capillary pressure is found to be negligible compared to the pressure developed within the bubble (see Fig. 4), even prior to the initial fracturing event (t = 444 s). It subsides even more with time as fracturing proceeds and the bubble grows.

Fig. 7. The evolution of crack increment, a, at the bubble front shown as a function of the vertical coordinate z. The time period from the first fracturing event (t = 444 s) until the moment when a at the tail of the crack drops to zero (t = 3524 s) is presented on the picture. The top curve corresponds to the initial (t = 444 s) distribution of the crack increment at the front when z A = −0.004 m, z B = 0.004 m (vertical coordinates of the bubble tail and head). The bottom curve corresponds to final moment of the fracturing, when a at the bubble tail drops to zero (t = 3524 s) and z A = −0.039 m, z B = 0.09 m. This figure demonstrates that crack propagation increment always remains constant at the bubble head and subsides at the bubble tail. After t = 3524 s bubble never grows at its tail.

section shows a closure of the bubble tail (w n = 0) thus indicating the onset of its ascent. It is demonstrated that a bubble with an initial penny-shaped configuration (R = 4 mm, 1.3 mm of the equivalent spherical radius), and with an almost elliptic opening profile, became highly asymmetrical as shown in the front view, with an opening resembling an “inverted tear drop”. The bubble equivalent spherical radius reaches 11 mm just prior to its rise. In order to study the bubble size and shape evolution, we track the ratio of the modeled bubble surface area to its volume (S / V ) with time (Fig. 9(a)). Furthermore, we calculate this parameter (S / V ) for the equivalent oblate-spheroidal bubble with

342

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

Fig. 8. The initial (t = 444 s) and final (t = 3524 s) bubble shape and size in the front view (panel a, X–Z plane in Fig. 1), and in the cross-section (panel b, Z –Y plane in Fig. 1). Both panels show that bubble grew significantly and became asymmetrical with respect to the initial penny-shaped configuration. The cross-section shows a closure of the bubble tail (w n = 0) thus indicating the onset of its ascent. The figure demonstrates that a bubble with an initial penny-shaped configuration (R = 4 mm, 1.3 mm of the equivalent spherical radius), and with an almost elliptic opening profile, became highly asymmetrical, with an opening resembling an “inverted tear drop”. The bubble equivalent spherical radius reaches 11 mm just prior to its rise.

Fig. 9. (a) Evolution of the bubble surface area to its volume ratio (S / V ) depicted for the modeled bubble, for the equivalent oblate-spheroidal bubble, and for the equivalent spherical bubble (see the text for the details). It is demonstrated that S / V ratio of the modeled and of the equivalent oblate-spheroidal bubbles remain very close and subside with time. S / V ratio of the equivalent spherical bubble also subsides with time, being much smaller than that of the modeled bubble. (b) Evolution of the aspect ratio (a/c, the relation of the major semi-axis to the minor one for the equivalent oblate-spheroidal bubble) grows with time.

Fig. 10. Allometric changes in the bubble proportions. For the isometric bubble growth pattern, the relation S / V ∼ S −1/3 holds during the entire growth period. However, for the modeled bubble a slope of S / V vs. S (depicted in the log scale) increases from about −0.22 to −0.04 with bubble growth. This demonstrates that the bubble grows allometrically, and its surface area remains relatively higher than expected from the isometric growth pattern. The changes in the bubble shape increase the effectiveness of the solute supply from pore water to the bubble interior.

be represented by the equivalent oblate-spheroidal bubble with no buoyancy affecting its shape. In this case the calculated equivalent bubble opening, c, compensates for the measure of the apparent bubble asymmetry. Just for the sake of comparison, the S / V ratio of the equivalent spherical bubble (with volume identical to the modeled one) also subsides with time as the bubble grows, being much smaller than S / V of the modeled bubble. In order to get an insight into the evolution of the bubble shape with its growth, one can evaluate allometric changes in its proportion (Thompson, 1992). If the bubble would grow isometrically with respect to its initial shape, its surface area would scale with its linear dimension as S ∼ l2 , and the volume as V ∼ l3 . As such, S / V ∼ S −1/3 during the entire period of bubble growth (Fig. 10). However, a slope of S / V vs. S for the modeled bubble is higher than −1/3 and increases from about −0.22 to −0.04 with bubble growth. This demonstrates that the bubble grows allometrically, and its surface area remains relatively higher than expected from the isometric growth pattern. The changes in the bubble shape certainly increase the effectiveness of the solute supply from pore water to the bubble interior as the bubble grows (occurred for the large bubble within a very short time period, Fig. 5). 4. Discussion 4.1. Bubble size and its controls

volume identical to the modeled one. This equivalent bubble is presented by a penny-shaped crack under no buoyancy (Fig. 1) with the radius equal to the half the difference in the modeled bubble head and tail coordinates (a = ( z B − z A )/2), z B = −z A at each time moment (e.g. Fig. 8). The minor bubble semiaxis (maximal bubble opening) is calculated as c = V /(4/3π a2 ) (Beyer, 1987), where V is modeled bubble volume at each moment. The surface area of this equivalent bubble is calculated as S = 2π a2 (1 + 1−ee tanh−1 (e )), e 2 = 1 − ac 2 (Beyer, 1987). Fig. 9(a) demonstrates that S / V ratio of the modeled and the equivalent oblate-spheroidal bubbles remain very close and subside with time, when a/c aspect ratio of the equivalent oblate-spheroidal bubble grows (Fig. 9(b)). Therefore, the asymmetric bubble can 2

2

Studying the bubbles within the sediment by using CT scans, Wilkens and Richardson (1998), Anderson et al. (1998), and Abegg and Anderson (1997) discerned that the equivalent spherical radius of the bubbles at the studied location and depth ranged mainly between 0.4 mm and 5 mm. However, individual bubbles with radius of ∼10 mm were also observed (Anderson et al., 1998). The size of a modeled “mature” bubble prior to its rise (∼11 mm of equivalent spherical radius) exceeded the main span of the observed bubble dimensions, approaching the largest observed bubbles size. This phenomenon is discussed below. The ambient concentration of the dissolved methane prescribed in this study according to the field data (Section 2.3) allows a permanent solute supply to the growing bubble. Hence, under

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

these ambient conditions, the bubble size and shape become entirely controlled by the mechanical properties of the sediment. This means that the bubble will grow to its maximal size avoiding no growth conditions caused by the insufficient solute supply to the bubble discussed in Algar and Boudreau (2010). Once the bubble terminal dimensions are reached (and when the SIF at its tip remains in close vicinity to the critical value), the bubble “cuts” the sediment by the process resembling a dynamic fracturing, and is eventually released from the sediment-water interface. This process may even be accelerated by the decrease in the critical SIF in the sub-bottom mud in reverse proportionality to the depth. This could be a reason why the large bubbles with dimensions similar to the modeled one are rarely observed (Wilkens and Richardson, 1998; Anderson et al., 1998; and Abegg and Anderson, 1997). However, inhomogeneities in material properties (especially in SIF) may stop the bubble growth and subsequent propagation or change their path. In this situation the large S / V ratio of the bubble and its increasing allometry (Figs. 9, 10) may play a crucial role by allowing the bubble dissolution close to the sediment-water interface due to the local methanotrophy. This may induce retardation of the bubble migration and even its partial dissolution. The viscous aspects of the sediment behavior are assumed to be negligible. Later we wish to discuss how the material parameters of the muddy sediment affect the terminal bubble size. The critical stress intensity factor, K Ic at depths 80–150 cm in muddy sediments of Eckernförde Bay was predicted to range between 600 and 1300 Pa m1/2 (Algar and Boudreau, 2010) and agreed with values estimated by Johnson et al. (2002). K Ic measured by Johnson et al. (2012) in Nova Scotia aquatic surface sediments by various methods spanned a wide range from zero (near the surface) to about 1800 Pa m1/2 (at a sediment depth of 0.2 m), changing inversely proportional to sediment porosity. In contrast, in gelatin, taken as a muddy sediment surrogate, K Ic varied between 28 to 300 Pa m1/2 being dependent on “porosity” (i.e., percent of water by mass; Menand and Tait, 2001; Rivalta and Dahm, 2006; Barry et al., 2010). To define a dependence of the bubble radius on critical SIF, Algar et al. (2011b) derived an analytical expression defining the critical radius of the buoyant penny-shaped crack at the beginning of its rise as proportional to the critical SIF to the power of 2/3. We verify this hypothesis by our numerical modeling and analytical derivations that account for the differential fracturing over the bubble front as it naturally occurs during the bubble growth. For a constant Young’s modulus and Poisson ratio, only the displacement component near the crack front normal to its surface, w nP , determines the value of K I distributed at the crack front (Eq. (8)). That is the difference in SIF between the crack head and any other point at the crack front, K I , stays constant as bubble inflates, which follows from the principles of linear elasticity. For instance, for the initial penny-shaped bubble under the buoyancy, this difference K IAB = 2 Pa m1/2 (Fig. 3) remains constant as SIF at the bubble head, K I B , grows during the initial elastic bubble inflation with no dependence on its critical value at the fracturing, K Ic . Further, assuming that K I B = K I max = K Ic at the crack head at each fracturing event and rearranging Eq. (9), yields:

a = amax (1 − K I / K Ic )2

(11)

where a is the local crack increment at the bubble front, and amax is maximal crack increment at the bubble head. Eq. (11) shows that a stays constant at the crack head ( K I = 0). When the fracturing starts at the earlier stage, i.e. when K Ic is smaller, a subsides faster toward the crack tail and becomes zero when

343

K I = K Ic . This results in the smaller terminal bubble size just prior to its rise. Therefore, to agree with the representative equivalent radii of the bubbles measured by Wilkens and Richardson (1998), Anderson et al. (1998), Abegg and Anderson (1997), the critical SIF must be notably smaller than 200 Pa m3/2 suggested in the literature and used in this study. In addition, in nature when the growth of bubbles occurs in sediments with a net of open or temporally closed ducts of low SIF (Chanton et al., 1989; Martens and Klump, 1980; Algar et al., 2011a), the bubbles will grow and rise at significantly lower terminal sizes than those obtained in our simulations. 4.2. Bubble shape The shape of the growing bubble (Fig. 8(a)) shows an asymmetry in the bubble front view with respect to its initial pennyshaped configuration. We found that this apparent asymmetry can also be observed in photos of the experiments on bubble grown in gelatin (Barry et al., 2010, Fig. 4(a)), and in soft muddy sediment (Boudreau et al., 2005, Fig. 2). This asymmetry has never been addressed in the recent studies. Our simulations that allow a differential fracturing over the bubble front succeeded to capture this feature. It follows from the above derivations (see Eq. (11)) that the measure of this apparent asymmetry depends on the critical SIF. The bubble allometry depending in turn on the measure of the bubble asymmetry increases the effectiveness of the solute supply to the growing bubble via its surface. Moreover, the bubble cross-section prior to its rise resembles an inverted tear drop (Fig. 8(b)) in agreement with earlier predictions of Weertman (1971), Nunn (1996), Nunn and Meulbroek (2002), Takada (1990), Algar et al. (2011b). 4.3. Predicting the delivery of gaseous methane to the atmosphere Bubbles are an important methane carrier to the atmosphere. In contrast to the dissolved methane, methane bubble can bypass processes of oxidation in upper sediment layers due to their fast rise velocity. A portion of the bubble methane is dissolved into the water column after the bubble release from the sediment. In the water, bubble rise velocity and gas transfer across the bubble surface are primarily affected by the bubble size (Clift et al., 1978; Leifer and Patro, 2002; McGinnis et al., 2006; Ostrovsky et al., 2008). Specifically, it was suggested that the bubble methane fraction transported to the atmosphere may be approximated by a power law of the equivalent spherical radius of the bubble emitted from the sediment (Leifer and Patro, 2002). In our case, the bubbles released from the sediment-water interface with equivalent spherical radius of r = 11 mm at the shallow water depth of 26 m would bring more than 90% of its initial methane to the atmosphere (Leifer and Patro, 2002; McGinnis et al., 2006). In contrast, bubbles with r = 5 mm would bring about 50% only (Ostrovsky et al., 2008). Therefore, the precise modeling of the bubble size and shape evolution in muddy sediments is of the major importance. 4.4. Fracturing of the fine-grained materials and other implications In various geological systems, including also the modeled one, the buoyancy-driven hydro-fracturing controls the fluid relative motion against solids. For instance, pervasive opening-mode subvertical fracturing observed in the low-permeability sandstone (The Upper Cretaceous Mesaverde Group in the Piceance Basin, Colorado) emerges under the high pore-fluid pressure in response to gas generation during the rock charging (Fall et al., 2012). On

344

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

a much larger scale, buoyancy-induced hydro-fracturing drives the magma intrusions and the dikes propagation (e.g. Lister, 1990). The fracture growth mechanism in the various fine-grained materials reveals a remarkable similarity. Macroscopic opening-mode fracture propagation direction controlled by the global stress field is demonstrated in the different systems: In the soft muddy sediment it is driven by gas pressure evolution (Boudreau, 2012), while in the diatomaceous mudstone it is induced by clinker and chert spheroids growth resulting in mineralogical and textural host rock reorganization (Eichhubl, 2004), etc. In the microscopic scale fracture growth can be explained by pore growth and coalescence. While the opening-mode fracturing occurs just in front of the macroscopic fracture tip, it was suggested by Eichhubl (2004) that in many substances (metals, ceramics, clinker and chert growth in the mudstone) creep-induced damage accumulates also obliquely to the main fracturing direction, driven by the shear stress concentrations. Therefore, the finding that on the microscopic scale the fracture should propagate in a zig-zag mode around the main macroscopic fracture propagation direction (Eichhubl, 2004) may also be applicable to the studied bubble growth in muddy sediment. An important question remaining is if the vertically-aligned fractures that could clearly be associated with bubbles-migrations paths in the muddy sediment and were numerously reproduced in the experiments (Boudreau, 2012 and references therein), have ever been identified in the mudrocks. Along with inter- and intra-granular cementation persisting in the mudrocks (Milliken and Day-Stirrat, 2013), the studies discern bedding-parallel planar deformation bands with a quantifiable compactional component, bioturbation-related particle alignments (Milliken and Reed, 2010), spherical poorly-cemented carbonate concretions (Raiswell and Fisher, 2000), and some other localizations. However, subvertical fractures in the mudrocks that could clearly be associated with bubbles-migrations paths, have never been discussed in the literature. We suggest that because these fractures are closed most of their lifetime by non-zero remote horizontal stresses, they experience a partial healing. Moreover, their surfaces that remain flat (zero curvature) also during the burial, are defined with chemical potential identical to that in the adjacent locations (Eichhubl, 2004; Lasaga, 1997). This avoids a preferential dissolution/precipitation on these fractures surfaces in contrast to what could be expected in the open fractures/veins/pores characterized by the non-zero curvature (Katsman, 2010). The observed preferential bubbles release under the low tides or water level (Chanton et al., 1989; Martens and Klump, 1980; Algar et al., 2011a) points onto the long-time fractures existence. Permeability evaluation of the fractured muddy sediment still remains scarce. However, modeling of other systems, e.g. of the fractured low-permeable carbonate reservoirs, suggests that even a poorly connected, widely spaced (∼10 m) system of the fractures increases the permeability two- to tenfold as compared with unfractured host material (e.g. Philip et al., 2005). Therefore, it is of the major importance for the wide interdisciplinary community to properly qualify methane bubbles growth and migration in the low-permeable muddy aquatic sediments. 5. Conclusions In this paper we present the first numerical modeling of the buoyancy-driven methane bubble growth in fine-grained muddy sediment prior to the beginning of its rise. 1. We designed a coupled mechanical/reaction-transport numerical model that enables a differential fracturing over the bubble front (as it occurs in nature). The model simulates a fracturing

2.

3.

4.

5.

6.

7.

8.

increment that stays constant at the bubble head and subsides towards bubble tail during its growth. We show that the differential fracturing over the bubble front controls the bubble shape and size temporal evolution. It is significantly affected by the critical stress intensity factor of the muddy sediment. Capillary pressure at the bubble surface is negligible even for the initial penny-shaped bubble with respect to the pressure developed within the bubble. Our simulations reveal a high asymmetry in the bubble shape growing with time, with respect to its initial symmetric penny-shaped configuration. We also confirmed the earlier predictions about the “inverted tear-drop” bubble cross-section just prior to the beginning of its rise. The duration of the elastic expansion stage and the value of the fracturing-induced pressure drop increment during the bubble growth subside with time as bubble approaches its terminal size (prior to its rise). We show that the allometry in the bubble amplifies with time, which affects the effectiveness of the solute supply to the growing bubble. We verify some hypotheses suggested by the previous studies and obtain the new insights based on our extensive simulations. Simulations of the terminal bubble dimensions allow predicting the delivery of gaseous methane from the sediment to the atmosphere via the water column.

Acknowledgements We wish to thank COMSOL Multiphysics support team for technical assistance, D.F. McGinnis for helpful comments, and three anonymous reviewers for constructive criticism that helped improving the clarity and quality of our presentation. This study was supported by grant from the Ministry of Energy and Water Resources, Israel (No. 29-17-033). References Abegg, F., Anderson, A.L., 1997. The acoustic turbid layer in muddy sediments of Eckernfoerde Bay, Western Baltic: Methane concentration, saturation and bubble characteristics. Mar. Geol. 137 (1–2), 137–147. Algar, C.K., Boudreau, B.P., 2009. Transient growth of an isolated bubble in muddy, fine-grained sediments. Geochim. Cosmochim. Acta 73, 2581–2591. Algar, C.K., Boudreau, B.P., 2010. Stability of bubbles in a linear elastic medium: Implications for bubble growth in Marine sediments. J. Geophys. Res. 115, F03012. Algar, C.K., Boudreau, B.P., Barry, M.A., 2011a. Release of multiple bubbles from cohesive sediments. Geophys. Res. Lett. 38, L08606. Algar, C.K., Boudreau, B.P., Barry, M.A., 2011b. Initial rise of bubbles in cohesive sediments by a process of viscoelastic fracture. J. Geophys. Res. 116, B04207. Anderson, A.L., Abegg, F., Hawkins, J.A., Duncan, M.E., Lyons, A.P., 1998. Bubble population and acoustic interaction with the gassy floor of Eckenforde Bay. Cont. Shelf Res. 18, 1807–1838. Barry, M.A., Boudreau, B.P., Johnson, B.D., Reed, A.H., 2010. First order description of the mechanical fracture behavior of fine grained surficial marine sediments during gas bubble growth. J. Geophys. Res. 115, F04029. Best, A.I., Tuffin, M.D.J., Dix, J.K., Bull, J.M., 2004. Tidal height and frequency dependence of acoustic velocity and attenuation in shallow gassy marine sediments. J. Geophys. Res. 109, B08101. Beyer, W.H., 1987. CRC Standard Mathematical Tables, 28th ed. CRC Press, Boca Raton. Boudreau, B., 1997. Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments. Springer-Verlag, Berlin, 414 pp. Boudreau, B.P., 2012. The physics of bubbles in surficial, soft, cohesive sediments. Mar. Petroleum Geol. 38, 1–18. Boudreau, B.P., Algar, C., Johnson, B.D., Croudace, I., Reed, A., Furukawa, Y., Dorgan, K.M., Jumars, P.A., Grade, A.S., Gardiner, B.S., 2005. Bubble growth and rise in soft sediments. Geology 33 (6), 517–520. Broek, D., 1986. Elementary Engineering Fracture Mechanics. Kluwer Academic, Dordrecht, Netherlands. Buckingham, M.J., 2005. Compressional and shear wave properties of marine sediments: Comparisons between theory and data. J. Acoust. Soc. Am. 117 (1), 137–152.

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

Chanton, J., Martens, C., Kelley, C.A., 1989. Gas transport from methane-saturated, tidal freshwater and wetland sediments. Limnol. Oceanogr. 34 (5), 807–819. Choi, J.-H., Seol, Y., Boswell, R., Juanes, R., 2011. X-ray computed-tomography imaging of gas migration in water-saturated sediments: From capillary invasion to conduit opening. Geophys. Res. Lett. 38, L17310. Citarella, R., Cricri, G., 2010. Comparison of DBEM and FEM crack path predictions in a notched shaft under torsion. Eng. Fract. Mech. 77, 1730–1749. Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops, and Particles. Academic Press, New York, 380 pp. Couroneau, N., Royer, J., 2000. Simplifying hypotheses for the fatigue growth analysis of surface cracks in round bars. Comput. Struct. 77 (4), 381–389. Dale, A.W., Aguilera, D.R., Regnier, P., Fossing, H., Knab, N.J., Jørgensen, B.B., 2008. Seasonal dynamics of the depth and rate of anaerobic oxidation of methane in Aarhus Bay (Denmark) sediments. J. Mar. Res. 66, 127–155. Davie, M.K., Buffett, B.A., 2001. A numerical model for the formation of gas hydrate below the seafloor. J. Geophys. Res. 106 (B1), 497–514. Davis, A.M., Huws, D.G., Haynes, R., 1996. Geophysical ground-truthing experiments in Eckernförde Bay. Geo Mar. Lett. 16, 160–166. Del Sontro, T., McGinnis, D.F., Sobek, S., Ostrovsky, I., Wehrli, B., 2010. Extreme methane emissions from a Swiss hydropower reservoir: Contribution from bubbling sediments. Environ. Sci. Technol. 44, 2419–2425. Dorgan, K.M., Arwade, S.R., Jumars, P.A., 2007. Burrowing in marine muds by crack propagation: kinematics and forces. J. Exp. Biol. 210, 4198–4212. Duan, Z., Moller, N., Greenberg, J., Weare, J.H., 1992. The prediction of methane solubility in natural waters to high ionic strength from 0 degrees to 250 degrees and from 200 to 1600 bar. Geochim. Cosmochim. Acta 56, 1451–1460. Eichhubl, P., 2004. Growth of ductile opening-mode fractures in geomaterials. In: Cosgrove, J.W., Engelder, T. (Eds.), The Initiation, Propagation, and Arrest of Joints and Other Fractures. In: Geological Society of London Special Publication, vol. 231, pp. 11–24. Fall, A., Eichhubl, P., Cumella, S.P., Bodnar, R.J., Laubach, S.E., Becker, S.P., 2012. Testing the basin-centered gas accumulation model using fluid inclusion observations: Southern Piceance Basin, Colorado. Am. Assoc. Pet. Geol. Bull. 96 (12), 2297–2318. Favelukis, M., 2004. Dynamics of foam growth: Bubble growth in a limited amount of liquid. Polym. Eng. Sci. 44 (10), 1900–1906. Gardiner, B.S., Boudreau, B.P., Johnson, B.D., 2003a. Growth of disk-shaped bubbles in sediments. Geochim. Cosmochim. Acta 67 (8), 1485–1494. Gardiner, B.S., Boudreau, B.P., Johnson, B.D., 2003b. Slow growth of an isolated diskshaped bubble of constant eccentricity in the presence of a distributed gas source. Appl. Math. Model. 27, 817–829. Gross, E., Seelig, T., 2006. Fracture Mechanics: With an Introduction to Micromechanics. Springer, Berlin. Guo, W., Shen, H., Li, H., 2003. Stress intensity factors for elliptical surface cracks in round bars with different stress concentration coefficient. Int. J. Fatigue 25, 733–741. Haeckel, M., Suess, E., Wallmann, K., Rickert, D., 2004. Rising methane gas bubbles form massive hydrate layers at the seafloor. Geochim. Cosmochim. Acta 68, 4335–4345. Haeckel, M.P., Boudreau, B., Wallmann, K., 2007. Bubble-induced porewater mixing: a 3-D model for deep porewater irrigation. Geochim. Cosmochim. Acta 71 (21), 5135–5154. Irwin, G.R., 1957. Analysis of stresses and strains near the end of a crack transversing a plate. J. Appl. Mech. 24, 361–364. Jain, A.K., Juanes, R., 2009. Preferential mode of gas invasion in sediments: Grainscale mechanistic model of coupled multiphase fluid flow and sediment mechanics. J. Geophys. Res., Solid Earth 114, B08101. Jiang, L., LeBlond, P.H., 1993. Numerical modeling of an underwater Bingham plastic mudslide and the waves which it generates. J. Geophys. Res. 98 (C6), 10303–10317. Johnson, B.D., Boudreau, B.P., Gardiner, B.S., Maass, R., 2002. Mechanical response of sediments to bubble growth. Mar. Geol. 187, 347–363. Johnson, B.D., Barry, M.A., Boudreau, B.P., Jumars, P.A., Dorgan, K.M., 2012. In situ tensile fracture toughness of surficial cohesive marine sediments. Geo Mar. Lett. 32, 39–48. Kassir, M.K., Sih, G.C., 1975. Three-dimensional Crack Problems. Mechanics of Fracture, vol. 2. Noordhoff international Publishing, Leyden. Katsman, R., 2010. Extensional veins induced by self-similar dissolution at stylolites: analytical modeling. Earth Planet. Sci. Lett. 299, 33–41. Keller, M., Stallard, R.F., 1994. Methane emission by bubbling from Gatun Lake, Panama. J. Geophys. Res. 99 (D4), 8307–8319. Lasaga, A.C., 1997. Kinetic Theory in the Earth Sciences. Princeton University Press, Princeton. NJ. Laubach, S.E., Eichhubl, P., Hilgers, C., Lander, R.H., 2010. Structural diagenesis. J. Struct. Geol. 32 (12), 1866–1872. Lawn, B.R., Wilshaw, T.R., 1975. Fracture of Brittle Solids. Cambridge University Press, New York. Leifer, I., Patro, R.K., 2002. The bubble mechanism for methane transport from the shallow sea bed to the surface: A review and sensitivity study. Cont. Shelf Res. 22, 2409–2428.

345

L’Esperance, J.C., Boudreau, B.P., Barry, M.A., Johnson, B.D., 2012. Small-scale, highprecision and high-accuracy determination of Poisson’s ratios in cohesive marine sediments. Geo Mar. Lett. 33, 75–81. Lister, J.R., 1990. Buoyancy-driven fluid fracture: the effects of material toughness and of low-viscosity precursors. J. Fluid Mech. 210, 263–280. Malvern, L.E., 1969. Introduction to the Mechanics of a Continuous Medium. Prentice Hall, Inc., Englewood Cliffs, NJ. Martens, C.S., Klump, J.W., 1980. Biogeochemical cycling in an organic rich coastal marine basin. I. Methane sediment-water exchange processes. Geochim. Cosmochim. Acta 44, 471–490. Martens, C.S., Albert, D.B., Alperin, M.J., 1998. Biogeochemical processes controlling methane in gassy coastal sediments – Part 1. A model coupling organic matter flux to gas production, oxidation and transport. Cont. Shelf Res. 18, 1741–1770. McGinnis, D.F., Greinert, J., Artemov, Y., Beaubien, S.E., Wuest, A., 2006. Fate of rising methane bubbles in stratified waters: How much methane reaches the atmosphere? J. Geophys. Res. 111, C09007. Mei, C.C., Liu, K.F., 1987. A Bingham-plastic model for a muddy seabed under long waves. J. Geophys. Res. 92 (C13), 14581–14594. Menand, T., Tait, S., 2001. A phenomenological model for precursor volcanic eruptions. Nature 411, 7. Milliken, K.L., Day-Stirrat, R.J., 2013. Cementation in mudrocks: Brief review with examples from Cratonic Basin Mudrocks. In: Chatellier, J.-Y., Jarvie, D.M. (Eds.), Critical Assessment of Shale Resource Plays. In: AAPG Mem., vol. 103, pp. 133–150. Milliken, K.L., Reed, R.M., 2010. Multiple causes of diagenetic fabric anisotropy in weakly consolidated mud, Nankai accretionary prism, IODP Expedition 316. J. Struct. Geol. 32, 1887–1898. Mogollon, J.M., L’heureux, I., Dale, A.W., Regnier, P., 2009. Methane gas-phase dynamics in marine sediments: A model study. Am. J. Sci. 309, 189–220. Mogollon, J.M., Dale, A.W., L’heureux, I., Regnier, P., 2011. Impact of seasonal temperature and pressure changes on methane gas production, dissolution, and transport in unfractured sediments. J. Geophys. Res. 116, G03031. Navon, O., Chekhmir, A., Lyakhovsky, V., 1998. Bubble growth in highly viscous melts: Theory, experiments, and autoexplosivity of dome lavas. Earth Planet. Sci. Lett. 160, 763–776. Nunn, J.A., 1996. Buoyancy-driven propagation of isolated fluid-filledfractures: Implications for fluid transport in Gulf of Mexico geopressured sediments. J. Geophys. Res. 101 (B2), 2963–2970. Nunn, J.A., Meulbroek, P., 2002. Kilometer-scale upward migration of hydrocarbons in geopressured sediments by buoyancy-driven propagation of methane-filled fractures. Am. Assoc. Pet. Geol. Bull. 86 (5), 907–918. Oldenburg, C.M., Lewicki, J.L., Dobeck, L., Spangler, L., 2010. Modeling Gas Transport in the Shallow Subsurface During the ZERT CO2 Release Test. Transp. Porous Media 82, 77–92. Ostrovsky, I., McGinnis, D.F., Lapidus, L., Eckert, W., 2008. Quantifying gas ebullition with echosound: The role of methane transport by bubbles in medium-sized lake. Limnol. Oceanogr., Methods 6, 105–118. Philip, Z.G., Jennings, J.W., Olson, J.E., Laubach, S.E., Holder, J., 2005. Modeling coupled fracture-matrix fluid flow in geomechanically simulated fracture networks. SPE Reserv. Eval. Eng. 8 (4), 300–309. Pinder, G.F., Gray, W.G., 2008. Essentials of Multiphase Flow and Transport in Porous Media. John Wiley & Sons, Inc., Hoboken, NJ. Proussevitch, A.A., Sahagian, D.L., 1996. Dynamics of coupled diffusive and compressive bubble growth in magmatic systems. J. Geophys. Res. 101, 17447–17455. Proussevitch, A.A., Sahagian, D.L., Anderson, A.T., 1993. Dynamics of diffusive bubble growth in magmas: Isothermal case. J. Geophys. Res. 98 (22), 283–307. Raiswell, R., Fisher, Q.J., 2000. Mudrock-hosted carbonate concretions: A review of growth mechanisms and their influence on chemical and isotopic composition. J. Geol. Soc. 157, 239–251. Richardson, M.D., Briggs, K.B., 1996. In situ and laboratory geoacoustic measurements in soft mud and hard-packed sand sediments: Implications for highfrequency acoustic propagation and scattering. Geo Mar. Lett. 16, 196–203. Rivalta, E., Dahm, T., 2006. Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface. Geophys. J. Int. 166, 1424–1439. Rowland, F.S., 1985. Methane and chlorocarbons in the Earth’s atmosphere. Orig. Life 15, 279–297. Rubin, H., Yaniv, S., Spiller, M., Köngeter, J., 2008. Parameters that control the cleanup of fractured permeable aquifers. J. Contam. Hydrol. 96, 128–149. Seevee, J.E., 2010. Effective porosity measurement in marine clay. J. Environ. Eng. 136 (7), 674–681. Sih, G.C., 1974. Strain energy density factor applied to mixed mode crack problems. Int. J. Fract. 10, 305–321. Silva, A.J., Brandes, H.G., 1998. Geotechnical properties and behavior of highporosity, organic-rich sediments in Eckernförde Bay, Germany. Cont. Shelf Res. 18, 1917–1938. Sneddon, I.N., 1946. The distribution of stress in the neighborhood of a crack in an elastic solid. Proc. R. Soc. Lond. A 187, 229–260. Soumis, N., Lucotte, M., Canuel, R., Weissenberger, S., Houel, S., Larose, C., Duchemin, E., 2005. Hydroelectric reservoirs as anthropogenic sources of greenhouse gases. In: Lehr, J.H., Keeley, J. (Eds.), Water Encyclopedia: Surface and Agricultural Water. Wiley–Interscience, Hoboken, NJ, pp. 203–210.

346

R. Katsman et al. / Earth and Planetary Science Letters 377–378 (2013) 336–346

St. Louis, V.L., Kelly, C.A., Duchemin, E., Rudd, J.W.M., Rosenberg, D.M., 2000. Reservoir surfaces as sources of greenhouse gases to the atmosphere: A global estimate. Bioscience 50, 766–775. Tabiei, A., Wu, J., 2003. Development of the DYNA3D simulation code with automated fracture procedure for brick elements. Int. J. Numer. Methods Eng. 57, 1979–2006. Takada, A., 1990. Experimental study on propagation of liquid-filled crack in gelatin: Shape and velocity in hydrostatic stress conditions. J. Geophys. Res. 95, 8471–8481. Thompson, D.W., 1992. On Growth and Form. Cambridge University Press, Cambridge. Turcotte, D.L., Shubert, G., 1982. Geodynamics Applications of Continuum Physics to Geological Problems. Wiley, New York.

van Kessel, T., van Kesteren, W.G.M., 2002. Gas production and transport in artificial sludge deposits. Waste Manag. 22, 19–28. Weertman, J., 1971. Theory of water-filled crevasses in glaciers applied to vertical magma transport beneath ocean ridges. J. Geophys. Res. 76 (5), 1171–1183. Wheeler, S.L., 1988. A conceptual model for soils containing large gas bubbles. Geotechnique 38, 389–397. Wilkens, R.H., Richardson, M.D., 1998. The influence of gas bubbles on sediment acoustic properties: In situ, laboratory, and theoretical results from Eckernförde Bay, Baltic sea. Cont. Shelf Res. 18, 1859–1892. Winterwerp, J.C., van Kesteren, W.G.M., 2004. Introduction to the Physics of Cohesive Sediments in the Marine Environment. Elsevier, Amsterdam.