Transient growth of an isolated bubble in muddy, fine-grained sediments

Available online at www.sciencedirect.com

Geochimica et Cosmochimica Acta 73 (2009) 2581–2591 www.elsevier.com/locate/gca

Transient growth of an isolated bubble in muddy, fine-grained sediments C.K. Algar *, B.P. Boudreau Department of Oceanography, Dalhousie University, 1355 Oxford St., Halifax, NS, Canada B3H 4J1 Received 17 September 2008; accepted in revised form 10 February 2009; available online 20 February 2009

Abstract Methane bubbles in fine-grained sediments have been shown to grow initially by elastic expansion and fracture. A previous growth model assumed quasi-steady state diffusion in which the methane porewater concentration quickly adjusted to changes in bubble geometry [Gardiner B. S, Boudreau B. P and Johnson B. D. (2003a) Growth of disk-shaped bubbles in sediments. Geochim. Cosmochim. Acta, 67 (8), 1485–1494]. Here, we present a finite-element model that solves the transient form of the reaction–diffusion equation, and the coupled linear elastic fracture mechanics (LEFM). In so doing we also employ a new theory for the post-fracture bubble sizes, based upon the full principles of LEFM. Our findings indicate that the quasi-steady state assumption is flawed due to violation of conservation of mass during fracture events. When the new model is applied to sediment conditions found at Cape Lookout Bight, NC, USA, it is found that bubbles grow somewhat faster than previously thought. A reference bubble of 0.5 cm3 will form in about 6 days, 2.5 days quicker than the old model predicted. Moreover, typical bubbles of 0.04 cm3 for this site can grow in as little as a day and a half. We examined the sensitively of the finite-element model to the various parameters in order to gain an understanding of how bubbles may behave under different sediment conditions. The influence of tides on bubble growth, through the process of rectified diffusion, was also examined and it was found that this had little influence upon growth. Ó 2009 Elsevier Ltd. All rights reserved.

1. INTRODUCTION The degradation of organic matter in sediments results in the production of volatile substances, one of which is methane. If production rates are high, this gas has the potential to saturate the porewater and form bubbles. Dissolution of gas hydrates can also lead to bubble formation. Gassy sediments occur throughout continental margins (Martens and Klump, 1980; Wever et al., 1998) and can have a significant impact on the geotechnical and geochemical properties of the sediments. The presence of bubbles can destabilize sediments, resulting in slope failure and has implications for the construction of any structures anchored in the sediments (Sills and Wheeler, 1992). The presence of gas in sediment is also an impediment to acoustic imaging of the *

Corresponding author. E-mail address: [email protected] (C.K. Algar).

0016-7037/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.gca.2009.02.008

ocean bottom; specifically the shapes and sizes of gas bubbles affects the attenuation and scattering of acoustic signals and must be taken into account when imaging gassy sediments (Anderson and Hampton, 1980a,b). Large bubbles may rise through the sediment creating channels and tubes. These tubes irrigate sediments, facilitating the exchange of solutes with the overlaying porewater (Haeckel et al., 2007) and can resuspend sediments (Klein, 2006). Ebullition can represent a significant methane source to the water column and eventually the atmosphere (Judd, 2004). Methane is a potent greenhouse gas, and it has been shown from ice core records that atmospheric CH4 concentrations are closely tied to temperature (Brook et al., 1996); therefore, attempts to develop a global methane budget and its potential effect on global warming should take into account sediments, including wetlands, as a methane source. Despite the existence of gas bubbles in marine sediments and their influence on the surrounding environment, a

2582

C.K. Algar, B.P. Boudreau / Geochimica et Cosmochimica Acta 73 (2009) 2581–2591

quantitative mechanistic understanding of their growth and migration through sediments is only now beginning to emerge. According to Anderson et al., 1998 sediment gas bubbles can be classified into three different types. Type 1 bubbles are small interstitial bubbles entirely contained within the pore spaces of the solid sediment framework. Type 2 bubbles are gas reservoirs in which the bubble is larger then the pore size, but the sediment structure remains unaltered leaving large areas of ‘‘dry” sediment. Type 3 bubbles displaces the sediment framework to form sediment free, gas filled voids, and it is this type that is most common in the fine-grained muddy sediments common to areas of high methane production. The importance of grain size in determining the behavior of gassy sediments has recently been demonstrated by Jain and Juanes (2008) in their investigations of free gas migration at the base of the hydrated stability zone. They suggest that in course-grained sediments, gas uniformly invades the pore spaces (Type 2 bubbles), while movement in fine-grained sediments is governed by a fracture-dominated regime. Indeed, it has become apparent in the last few years that the Type 3 bubbles, common to muddy sediments, initially grow by way of elastic expansion and fracture of the surrounding sediments. This has been confirmed through observations using CT-scanning technology to view bubbles in sediment cores (Abegg and Anderson, 1997; Anderson et al., 1998; Boudreau et al., 2005), bubble growth experiments (Johnson et al., 2002) and modeling studies (Gardiner et al., 2003a,b). Bubbles that grow subsequently in locations where an initial bubble grew are governed by a complex set of dynamics not explored in this paper. Johnson et al. (2002) provided a framework in which to develop a mechanistic model of initial bubble growth, based upon the principles of linear elastic fracture mechanics (LEFM), whereby the fracture behavior is characterized by the fracture toughness parameter K 1c . Gardiner et al. (2003a) combined the principles of LEFM with a reaction–diffusion model to describe gas supply to the bubble and made estimates of bubble growth rates, predicting that bubbles grow to the large observed sizes in nature ð 0:5 cm3 Þ on a weekly timescale. A key assumption involved in the Gardiner et al. (2003a) model is that bubbles grow sufficiently slowly so that the solute concentration field next to the bubble can readily adjust to changes in bubble geometry, maintaining a quasisteady state. The reaction–diffusion portion of the model then reduces to finding the solution of the steady state diffusion equation in the presence of a constant distributed source, Dr2 w þ S ¼ 0;

ð1Þ

where D is the diffusion coefficient, w is the porewater methane concentration, S is the rate of methanogenesis and r2 is the Laplacian operator. This assumption of quasi-steady state allowed Gardiner et al. (2003b) to find an analytical solution to Eq. (1) in the oblate-spheroidal coordinate system by the method of separation of variables, in a manner similar to Lebedev (1972). This assumption, however, while central to the model, could not be explicitly tested.

Bubble growth, by its very nature is a transient process. During elastic expansion the concentration at the bubble surface continuously increases as gas diffuses into the bubble, and fracture events result in sudden increases in bubble volume and drops in internal gas pressure/concentration (Johnson et al., 2002). Both processes cause changes in the porewater concentration field, implying that the porewaters can never truly be at steady state and bringing into question the quasi-steady state assumption. If this assumption is not valid then the supply of gas to the bubble and, therefore, the growth rate would be different than predicted by Gardiner et al. (2003a). For this reason we present a finite-element model for LEFM/reaction–diffusion initial growth of bubbles that does not rely on a condition of quasi-steady state, but solves the true transient form, @w ¼ Dr2 w þ S: @t

ð2Þ

This allows us to test the validity of this assumption and obtain more accurate estimates of bubble growth in sediments. Before presenting this new solution we take a short departure in order to correct an error in the original presentation of the LEFM theory of bubble growth that has been carried through in the analysis of Gardiner et al. (2003a). 2. CRACK OPENING DISPLACEMENT In the LEFM theory, a bubble is treated as a pennyshaped crack imbedded in an infinite elastic medium. Growth occurs in two different phases. When the internal pressure is below a critical value, the bubble expands elastically by compressing and/or deforming the surrounding sediments. This does not change the length of the bubble’s long axis, a, but increases the width or short axis of the bubble, b (Fig. 1). When the internal bubble pressure exceeds a critical value (characterized by K 1c ), the bubble fractures the surrounding sediments, increasing the bubble length, and the width of the bubble then re-adjusts elastically. The elastic adjustment and expansion of the bubble walls is given by the crack opening displacement formula,

a COD

b

x

Fig. 1. A cross-section through an idealized oblate-spheroidal sediment bubble showing both the major axis ðaÞ, and minor axis ðbÞ. The distance separating the bubble walls is the crack opening displacement (COD) and is shown for one value of x, where x is the distance along the major axis. The minor axis is half the maximum crack opening displacement ðb ¼ COD=2Þ.

Bubble Growth in Sediments

COD ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4ðP b  P 0 Þ a2  x2 ; E

ð3Þ

where P b is the internal bubble pressure, P 0 is the ambient pressure, x is the distance along the major axis of the bubble and E is Young’s modulus. Eq. (3) is derived from fundamental elasticity theory for the case of an elliptical crack under conditions of plane stress (Broek, 1982). Plane stress conditions apply to geometries with plane dimensions much larger than their thickness and with forces applied only in the direction of the plane. However, this is not the loading condition experienced by a sediment bubble. A bubble experiences axially symmetric loading normal to its surface, with stresses and strains in all three dimensions. Sneddon (1946) gives the COD formula for the case of a penny-shaped void in a elastic solid as, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8ð1  m2 ÞðP b  P 0 Þ a2  x2 COD ¼ : ð4Þ pE This gives crack opening displacements that are roughly a factor of 2=p lower than the plane stress COD formula, resulting in the formation of bubbles which have lower aspect ratios ðb=aÞ than were previously predicted. Eq. (4) also introduces another parameter, Poisson’s ratio ðmÞ, into the model. Poisson’s ratio is the ratio between strain in one direction and strain in the other perpendicular directions. It normally ranges in value from 0 for a completely compressible material to 0.5 for an incompressible material. To date values of Poisson’s ratio in sediments are not well known. Although we would expect this value to approach that of an incompressible fluid for a fully saturated sediment it could decrease as the gas content of the sediment increases. The change in COD formula also changes a number of formulas in the LEFM bubble growth model. The expression for the critical pressure at which fracture occurs ðP c Þ, given by Eq. (16) in Johnson et al. (2002) becomes, Pc ¼

3=5 K 6=5 ð1  m2 Þ 1c p

121=5 ðEV b Þ1=5

:

2583

processes is that bubbles grow faster than previously thought: a reference bubble of 0.5 cm3 will form in about 7.7 days, a day and a half faster than predicted in Gardiner et al. (2003a) (Fig. 2). 3. FINITE-ELEMENT MODEL OF BUBBLE GROWTH With bubble mechanics now more consistent with LEFM theory, we turn our attention to the question of quasi-steady state. In the fully transient case (Eq. 2), a time dependency in the boundary condition at the bubble surface is introduced making this boundary condition inhomogeneous and the analytical solution becomes much more difficult; thus we implement a numerical method to find a solution. The finite-element method (FEM) is an ideal framework in which to approach this problem because of its ability to handle both curved boundaries and the two-way coupling between the reaction–diffusion and mechanics portions of the model. Furthermore the use of a numerical method eliminates the reliance on the rotational symmetry of a single bubble, and although we still exploit this symmetry in the results presented here, this need not be the case. In future work, we now have the ability to examine more complicated bubble orientations and boundary conditions that create porewater methane distributions more reflective of the natural environment. The model was developed and run in the commercially available finite-element modeling software COMSOL Multiphysics. What follows is a brief outline of our approach. The geometry for our model is the same as was employed by Gardiner et al. (2003a). A single isolated bubble of an oblate-spheroidal shape was embedded in a sphere of sediment with a radius of R1 . Due to axial symmetry only

ð5Þ

The propagation of a during a fracture event also changes from Eq. (15) in Johnson et al. (2002) to  1=3 3EðV b Þpf apf ¼ ð6Þ 16ð1  m2 ÞðP c  P 0 Þ where the subscript pf indicates post-fracture values. The net result of these changes is that bubbles will fracture at lower internal bubble pressure and that fractures, once initiated, will propagate further then previously thought. The question now becomes what affect these changes will have on the bubble growth times predicted by LEFM theory. The change in COD formula has two opposing effects upon the predicted growth rates. Eq. (4) shows that there will be less elastic expansion of the bubble walls for a given internal bubble pressure, favoring a slower growth rate. However, bubbles will also be more eccentric, giving a higher surface area to volume ratio, thereby increasing the diffusive flux of gas into the bubble, which favors a faster growth rate. The net effect of these competing

Fig. 2. Growth of a 0.5 cm3 reference bubble according to three different LEFM growth models; (1) the quasi-steady state model as is described in Gardiner et al. (2003a), (2) the quasi-steady state solution with the corrected COD formula (i.e. Eq. 4), and (3) the transient FEM model of LEFM growth described in Section 3. The model parameters for all three cases are those applied to Cape Lookout Bight (CLB) and are shown in Table 1. It can be seen that both the COD correction and the transient model decrease growth times compared to the original LEFM model.

2584

C.K. Algar, B.P. Boudreau / Geochimica et Cosmochimica Acta 73 (2009) 2581–2591

q 0.03

ð10Þ

where q is the bulk sediment density, c is the elasticity tensor as described in Timoshenko and Goodier (1970), and u is a vector giving the displacements in the radial and axial directions. For these simulations the effect of gravity is neglected and the internal bubble pressure above ambient is a load normal to the bubble surface,

R

1

Axial Coordinate (z)

@2u ¼ r  ðcruÞ @t2

n  cru ¼ wb RT  P 0 Sediment

ð11Þ

where R is the ideal gas constant and T is temperature. Finally V b is given by the following integral, Z Vb ¼ n  u dS; ðV bðt¼0Þ ¼ V 0 Þ: ð12Þ S¼a0

α

b

Bubble 0

Radial coordinate (r)

0.03

Fig. 3. Geometry for the FEM transient bubble growth model. The model uses the cylindrical (axisymmetric) coordinate system, specified by the coordinates r; z; / and represents 1 quarter of a slice through the bubble and surrounding sediment sphere. Due to rotational symmetry about the z axis the coordinate / does not come into the solution. ab is the bubble surface and R1 the radius of the surrounding sphere of sediment. Note that the size and proportions of the bubble have been exaggerate to make viewing easier.

one quarter of a single slice through the bubble and surrounding sediments needs to be modeled (Fig. 3). During elastic expansion the reaction–diffusion portion of the model calculates the porewater methane concentrations surrounding the growing bubble. This is done by solving Eq. (2), subject to the following boundary conditions: (1) at distances far from the bubble surface, the methane concentration approaches a constant value, wðR1 Þ ¼ w0

ð7Þ

where w0 is the ambient methane concentration and R1 is the radius of the model domain, and (2) at the bubble surface the porewater gas concentration is in equilibrium with the gas concentration inside the bubble according to Henry’s Law, wðab Þ ¼

wb ðtÞ k

where V 0 is the volume at the start of the elastic growth phase. According to Eq. (5), elastic growth continues until the overpressure reaches P c and fracture occurs. Fracture is treated as an instantaneous process, whereby the length of the crack ðaÞ extends according to Eq. (6) and the bubble walls re-adjust (Eq. 4). The methane porewater distribution at fracture is then mapped to its new post-fracture position and becomes the initial condition for the next phase of elastic expansion. 4. MODEL RESULTS In order to compare the results of the transient model to the quasi-steady state case, we performed simulations with the same parameter values used by Gardiner et al. (2003a). The parameters were chosen so as to mimic the Cape Lookout Bight (CLB) site in North Carolina (Table 1), where methanogenic rates have been measured (Crill and Martens, 1986); a more complete description of the site can be found in Boudreau et al. (2001). Our results indicate that the transient model predicts significantly faster growth than the quasi-steady state model (Fig. 2). Transient bubbles grow to a reference volume of 0.5 cm3 in just under 6.5 days, about 1 day faster than for the quasi-steady state case with the corrected COD formula (Eq. 4) and 2.5 days faster than in Gardiner et al. (2003a), an increase in growth rate of about 40%. The reason for the difference in growth rate is that, in the quasi-steady state model, the conservation of mass in the

ð8Þ

where ab represents the bubble surface, k is the dimensionless form of Henry’s Law constant and wb ðtÞ is the internal gas concentration of the bubble, which is a function of time. In addition, conservation of mass for the growing bubble can be stated as, Z @w @V b V b b þ wb n  Drw dS ð9Þ ¼ @t @t S¼ab where n is the outward normal at the bubble surface from the point of view of the sediments. During the elastic phase, the bubble volume ðV b Þ is determined from movement of the bubble walls as a result of the deformation of surrounding sediments according to the elasticity equation,

Table 1 Parameter values corresponding to conditions at Cape Lookout Bight (CLB), NC. Name

Symbol

Value

Diffusion coefficient Source strength Young’s modulus Poisson’s ratio Fracture toughness Temperature Ambient pressure Ambient methane Sediment radius

D S E m K 1c T P0 w0 R1

5  1010 m2 s1 3:47  106 mM s1 1:4  105 Nm2 0.3 300 Nm3=2 293 K 170; 000 Nm2 3.41 mM 0.03 m

Bubble Growth in Sediments

sediments surrounding a bubble during a fracture event is violated. This is shown in Fig. 4, which depicts the methane porewater concentration field both at the moment of fracture (a) and immediately after fracture (b) as calculated by the quasi-steady state model. During a fracture event, the amount of methane in a bubble remains constant but concentration of gas inside the bubble decreases due to an increase in volume. In quasi-steady state, this becomes the boundary condition for the initial steady state solution during the next phase of elastic growth. However, from Fig. 4 it is evident that under this assumption the sediments contain less methane post-fracture than was present immediately before fracture; mass has been lost due to the incorrect dynamics. The quasi-steady state model is based on the assumption that when the timescale of diffusion is much faster than the timescale of growth, the concentration field has time to adjust to changes in the bubble geometry. This assumption is flawed when applied to fracture events. The only way to produce the methane porewater distribution in Fig. 4b from a, is for the missing methane to flow across the outer boundary. This is physically impossible since it would require gas to move against the concentration gradient. In reality, there is an initial discontinuity at the surface of the bubble immediately after a fracture event. The concentration of gas inside the bubble decreases, but the porewater concentration at the bubble surface remains at its pre-fracture level. The ‘‘extra” gas in Fig. 4a compared to b diffuses into the growing bubble as the discontinuity disappears, but goes unaccounted for in the quasi-steady state model. This behavior is taken into account in the transient model, as is shown by the initially high flux in Fig. 5. It is this high flux across the bubble surface during the early stages of elastic expansion that accounts for the difference between the transient and quasi-steady state growth. These simulations clearly show that this error in conservation of mass is significant and cannot be ignored, as has previously been hypoth-

2585

Fig. 5. The rate of gas supply (mass transfer) to a growing bubble during the early stages of an elastic expansion phase for both the transient and quasi-steady state models. During the initial stages of an elastic expansion phase the flux is higher in the transient model (dashed line) then in the quasi-steady sate model (solid line) due to the initial discontinuity at the bubble surface. This accounts for the difference in growth rate between the two models At longer time periods the rate of gas supply for both cases becomes equal.

esized. The transient nature of the concentration field surrounding a growing bubble must therefore be taken into account. 5. INFLUENCE OF MODEL PARAMETERS This section examines the sensitivity of the model to the key parameters involved in bubble growth. In many cases, the variability and range of values of these parameters in natural sediments is not yet known; however, by examining a range of realistic parameter values we can gain

Fig. 4. Porewater methane concentration field surrounding a small ða ¼ 1 mmÞ bubble located at the origin for two stages of the growth cycle (a) pre-fracture and (b) post fracture, calculated using the quasi-steady state bubble growth model. The CLB (Table 1) set of parameters were used. The bubble is located in the lower left hand corner. It is clear that the methane concentration near the bubble is higher in the pre-fracture case compared to post-fracture, indicating violation of mass conservation in the quasi-steady state model.

2586

C.K. Algar, B.P. Boudreau / Geochimica et Cosmochimica Acta 73 (2009) 2581–2591

understanding into how bubbles grow under possible sediment conditions that may exist. Fig. 6 shows the effect of the source ðSÞ, Young’s modulus ðEÞ, poisson’s ratio ðmÞ, and fracture toughness ðK 1c Þ on growth rate. Values were chosen so that, to the best of our knowledge, they likely reflect the range of natural values. Methane production rates have been seen to vary over several orders of magnitude from around 1  108 mM s1 in low depositional environments (Albert et al., 1998) to upwards of 3  4  106 mM s1 for sediments with large inputs of organic matter, such as Cape Lookout Bight (Crill and Martens, 1986) or White Oak Estuary (Martens and Goldhaber, 1978; Chanton et al., 1989). Fig. 6a shows the effect of varying the source strength, S, and it clearly demonstrates the importance this parameter plays in determining how fast bubbles will grow. The greater the rate of methane production, the greater the rate of gas supply to the bubble, and the faster it will grow. A change in the model brought about by the use of the proper COD formula is the introduction of Poisson’s ratio as a parameter influencing growth mechanics. Finite-strain values of Poisson’s ratio in sediments are currently not well known. Fortunately this lack of knowledge is not a major problem since Poisson’s ratio has a secondary effect upon growth times (Fig. 6b). A larger Poisson’s ratio does increase the growth rate slightly, but this effect is small compared to the influence of the other mechanical parameters. The critical stress intensity factor, K 1c , is inversely related to growth rate as illustrated in Fig. 6c. Larger values of K 1c result in longer times between fracture events, thereby slowing down growth rates. During the early stage of an elastic expansion phase the flux of gas into the bubble is

high (Fig. 5), but slows down as the concentration gradient at the bubble surface is smoothed out. A low K 1c value results in a bubble that fractures early while the supply of gas to the bubble is still high, increasing the growth rate over a bubble that has a longer phase of elastic expansion (higher K 1c ). Young’s modulus, E, has an effect similar to that seen in the COD formula. A high Young’s modulus limits the amount of elastic expansion at a given overpressure, but at the same time results in bubbles with higher surface to volume ratios, allowing more surface for gas exchange. Simulations show that once again the surface area effect dominates (Fig. 6d), and bubbles with a higher Young’s modulus grow faster. These simulations show that a reference bubble of 0.5 cm3 will form on a timescale of a few days up to a week. This analysis indicates that the methane production rate is the key factor in determining how fast a bubble will grow, by ultimately controlling the supply of gas to the bubble. The mechanical properties are of less importance, but nevertheless cannot be ignored, since they are the main factor’s controlling bubble shape and thus the surface area over which gas may diffusive (see Gardiner et al. (2003a) for a good discussion on the comparison between oblate and spherical bubbles). Of these the fracture toughness ðK 1c Þ is of the most importance, with bubbles in weaker sediments growing faster than in stronger (higher K 1c ) sediments. Poisson’s ratio is the least important of these parameters; simulations show that the time required to grow a bubble of 0.5 cm3 only increases by around 10% over the entire range of possible values, 0 < m < 0:5.

Fig. 6. The influence of model parameters on bubble growth. Simulations were carried out over the expected range of natural variability for each parameter. In each plot, one of the parameters was varied while all others were held at the values in Table 1. The parameter varied in each case was; (a) source strength ðSÞ, (b) Poisson’s ratio ðmÞ, (c) fracture toughness ðK 1c Þ, (d) Young’s modulus ðEÞ. The discontinuities in the slopes of the growth curves clearly evident in (c) are real, and are, although less evident, present in all the growth curves. They are the result of the sudden increases in volume that occurs during fracture events when the crack length increases instantaneously.

Bubble Growth in Sediments

The FEM model of bubble growth conserves mass, but the question remains whether or not the model is consistent with observations of real sediment bubbles. Unfortunately such estimates, despite their importance in quantifying methane flux rates from sediments, remain rare. Growth times for individual bubbles in the natural environment have never been measured. Martens and Klump (1980) have reported bubbling rates at the sediment water interface of 220–1330 cm3 m2 h1 during low tidal periods at Cape Lookout Bight and the size distributions of bubbles escaping bubbling tubes (Table 2). However, these results represent a depth integrated result of an entire bubble population and take into account the mechanics of bubble rise, something our model does not yet include. Nevertheless it does provide us with a crude estimate of the timescale of bubble growth since we would expect that the bubbles being released must be forming on the order of a tidal cycle (12h) if the steady state methane porewater concentrations are to be maintained. Table 2 shows the growth times predicted by the LEFM bubble growth model for each bubble size. The parameters used are those listed in Table 1, where the source strength (S) of 3:47  106 mM s1 comes from the results of Crill and Martens (1986) who measured methane production rates at Cape Lookout Bight. Thirty-eight percent of the rising bubbles form in less than one complete tidal cycle, while the larger bubbles can take anywhere from one to several days to form. These estimates should however be viewed as maximum growth times, since the large bubbles are probably forming partly through movement and coalescence of smaller bubbles. Also the values of E and K 1c used are estimates of bulk sediment properties, and as Martens and Klump (1980) pointed out bubbles tend to rise from preexisting bubble tubes, inline with evidence that a fracture once formed is easier to reopen in the future (Boudreau et al., 2005). The mechanical properties of the bubble tubes are, therefore, probably different than those of the bulk sediment and favor faster growth rates. The LEFM-based bubble growth model predicts growth rates that seem to agree with those suggested for natural bubbles. Our previous non-resistant spherical model (Boudreau et al., 2001) is certainly now found to be severely wanting.

Table 2 Size of bubbles exiting bubble tubes during low tide bubbling periods as reported by Martens and Klump (1980) and growth times predicted by the transient LEFM growth model and the quasi-steady state model. Bubble radius range at 7.5 m depth (mm)

0.5–1 1–2 2–3 3–4 4–5 * **

Percentage of total bubbles

6 32 39 13 10

Mean volume per bubble (mm3)

1 7 38 101 210

A = transient finite-element LEFM model. B = quasi-steady state model.

Growth times A* (h)

B** ( h)

3.8 11 30 55 88

4.1 14 38 67 106

2587

Further evidence in support of our LEFM model comes from Klein (2006), who has carried out a series of laboratory experiments on the effect of bubble ebullition on sediment porewater exchange. In his experiments, bubbles were injected artificially into sediment cores resembling natural conditions at Barr Lake, Colorado, which is a seasonally anoxic eutrophic lake with corresponding methane production. Although bubbles were artificially injected, Klein reported the appearance of natural bubbles in the sediment cores appearing every two to three days, prior to the beginning of his experiment. Klein did not include any rates of methanogenesis or mechanical properties of his sediments; nevertheless, this does suggest a time scale for the formation of sediment bubbles consistent with our model predictions. 6. BUBBLE GROWTH AND TIDES Tidal cycles could potentially increase the rate of bubble growth through the process of rectified diffusion. Rectified diffusion describes the growth of a bubble due to an oscillating pressure field (Hsieh and Plesset, 1961). Such a pressure field causes a bubble to expand and contract. During contraction the gas concentration in the bubble rises above the ambient value and gas diffuses out of the bubble, while during expansion the bubble gas concentration decreases and gas diffuses into the bubble. This transfer is not balanced, but rather there is a net transfer of gas into the bubble due to the difference in surface areas during each half-cycle. For a spherical bubble in a non-resistant medium the mean rate of gas inflow due to the pressure oscillations is given by Hsieh and Plesset (1961) to be,  2   dm 8p MP ð13Þ ¼ DC 0 R dt 3 P0 were h::i denotes a time averaged quantity, MP is the amplitude of the pressure oscillations and R is the mean radius of the bubble. Using this, Boudreau et al. (2001) concluded that for spherical non-resistant bubbles, rectified diffusion due to tidal action does not contribute significantly to growth. The mechanical properties of the sediments cause the aspect ratio of a bubble to fluctuate over time due to elastic expansion and fracture events, complicating the application of rectified diffusion theory. To our knowledge such a situation has not been explored previously and no closed form solution exits. Fortunately the finite-element LEFM model can account for tidal effects by substituting the load on the bubble surface with a periodic boundary condition, n  cru ¼ wb RT  ðP 0 þ qw gA sin wtÞ

ð14Þ

where A is the tidal amplitude, w the period and qw the seawater density. In this model, the movement of the bubble wall is modeled explicitly and fully coupled to the reaction–diffusion portion of the model, and as a result, the rectified diffusion effect appears naturally. Martens and Klump (1980) stated that their study site at Cape Lookout Bight experiences an average tidal hight of 1 m in a water depth of 8 m, corresponding to ðMP =P 0 Þ2  0:05, indicating that

2588

C.K. Algar, B.P. Boudreau / Geochimica et Cosmochimica Acta 73 (2009) 2581–2591

the influence of rectified diffusion should be small. Fig. 7, which shows the accumulation of mass in a bubble during a model run with Cape Lookout Bight conditions, confirms this. However, the subtle influence of rectified diffusion upon growth is evident, with a greater increase in gas accumulation during low tides, when the bubble is expanding, than during high tides, when it is contracting. If these fluctuations did not result in a net influx of gas, then the no-tide growth curve (dashed line) would pass through the mean of the tidal bubble growth curve (solid line). This is not the case, but rather the tidal growth curve is shifted slightly higher, indicating a small increase in growth due to rectified diffusion. Fig. 8 compares this to what would be expected for a spherical bubble of similar volume according to Eq. (13) and shows that, while in both cases the influence of rectified diffusion is of similar order of magnitude, it is magnified in the LEFM case (by roughly a factor of 7). This is expected due to the increased surface area to volume ratio of an oblate spheroid compared to a sphere. A comparison with Fig. 7 shows that the accumulation of gas by rectified diffusion in the bubble is two orders of magnitude lower than gas accumulation due to all other process combined; consequently like the spherical case in Boudreau et al. (2001), the tidal influence on bubble growth is not significant. It is, however, important to point out that this result applies only to the initial growth of bubbles and says nothing about the influence of tides on the release and rise of bubbles into the water column. Many authors have reported ebullition occurring on tidal cycles (Martens and Klump, 1980; Jackson et al., 1998), but a release mechanism consistent with the mechanical response of the sediments has not yet been proposed. We plan to examine this in a future work.

Fig. 7. Accumulation of mass in a bubble vs. time during model runs both with and without the influence of tides. The tidal amplitude was taken to be 1 m as described in Martens and Klump (1980) and all parameters were set to CLB values (Table 1).

Fig. 8. Average net mass transfer due to rectified diffusion for an LEFM bubble growing under Cape Lookout Bight conditions (solid line) and compared to the case of spherical bubble of equal volume as predicted by Eq. (13) (dashed line). The effect of rectified diffusion in the LEFM case is greater then in the spherical case due to the increased surface area to volume ratio of an oblatespheroidal bubble.

7. FRACTURE PROPAGATION A current problem that still exists in the LEFM theory of bubble growth is that it says nothing about how far a fracture will propagate, only that it will propagate. Gardiner et al. (2003a) made note of this and in the absence of any set theory choose to treat each fracture event as an instantaneous pressure drop to a set value slightly above ambient ðP pf ¼ 1:01P 0 Þ. For consistency in comparing the two models we have used the same treatment in the simulations discussed so far. However, this fracture propagation rule creates a problem as a bubble grows to a large volume. Since, according to Eq. (6), P c is a function of bubble volume, as the bubble grows P c decreases and will eventually drop below P pf and the model breaks down. Furthermore, the decision to use a constant post-fracture pressure was one made for convenience and necessity, but it has been confirmed neither experimentally nor theoretically. It is for this reason that we propose a new criteria for fracture propagation based directly upon the principles of LEFM. To do this, we must first define the process zone. In LEFM theory, the process zone is the necessarily small region ahead of a crack, where the material deforms plastically and the elasticity equations no longer describe the stress field; rather the stress ahead of the crack tip is cut off at the yield stress ðry Þ, as is shown in Fig. 9. The plastic deformation of the process zone is what prevents the stress singularity that would otherwise be predicted by LEFM and is the region where the physical act of bond breaking takes place (Heimpel and Olson, 1994). We will now assume that upon fracture the material located inside the process zone fails instantaneously and the crack propagates the length of this zone, after which the internal pressure drops and growth is arrested. Such an approach has been used

Bubble Growth in Sediments

Fig. 9. The process zone and stress distribution surrounding a crack. The processes zone is a circle of material surrounding the crack tip where permanent plastic deformations take place and the elasticity equations no longer describe the stress distribution in the surrounding sediments. This zone is small compared to the crack length, but has been exaggerated here for clarity. The solid line depicts the stress distribution in the sediments along the crack plane as the crack tip is approached. The dashed line shows the stress distribution that would be predicted by LEFM theory inside the process zone if plastic deformations did not occur.

when applying fracture mechanics to the fracture of rocks during dyke formation (Rubin, 1993; Heimpel and Olson, 1994). There are many different formula for calculating the size and shape of the process zone, but one of the simplest and most widely accepted is the Dugdale zone (Gross and Seelig, 2006), which states that the plastic zone thickness ðDaÞ is given by,  2 p K1 Da ¼ ð15Þ 8 ry where ry is the yield strength of the material and K 1 is the stress intensity factor. By rearranging and substituting K 1c for K 1 we arrive at an expression giving the increase in crack length during a fracture event. This equation shows that during fracture the crack always increases by the same amount, regardless of the size of the bubble, and is governed by the ratio of the critical pressure and the yield strength of the sediment ðK 1c =ry Þ.

2589

This result underlines the necessity of determining the yield strength of soft sediments, something that is rarely, if ever, done. In the absence of such numbers, we will once again resort to modeling to understand the potential influence this parameter will have on the growth of bubbles. As a general rule, the theoretical yield strength of a solid with no flaws is approximately E=10 (Lawn and Wilshaw, 1975). Whereby a flaw is any microcrack or imperfection in the material which can serve to concentrate stress. This provides an upper estimate for the yield strength of the sediments. In reality, the yield strength is probably much lower since sediments, most certainly, contain many flaws. As well, their lack of tensile strength compared to more typical solids suggests that this ratio should be lower. To understand the extent to which fracture propagation and the new parameter, ry , influences bubble growth model runs were carried out with ry values at the theoretical E=10 ratio and 2/3 and 1/3 this value, which for the CLB site corresponds to ry ¼ 1:5  104 ; 1:0  104 and 5:0  103 Nm2 , respectively. Fig. 10 demonstrates that the internal bubble pressure can show different behavior depending upon the length of the process zone. Fig. 10a depicts a situation of high yield strength, ðry ¼ 1:5  104 Nm2 Þ. The result is a decreasing pressure record that appears to approach a flat-line trend. A smaller yield strength value ðry ¼ 5:0  103 Nm2 Þ produces a pressure record which is similar to the constant post-fracture pressure results reported in Gardiner et al. (2003a), whereby large discreet pressure drops clearly indicate fracture events. Interestingly both types of records are common during bubble growth experiments performed in sediments (Johnson et al., 2002). Previously flat-line pressure records in sediments were a bit of a mystery and we thought that these bubbles were growing by a different mechanism; our results clearly now show that both types of behavior are entirely consistent with LEFM theory. Furthermore considering the heterogeneous nature of sediments and the wide

Fig. 10. Internal pressure of a growing sediment bubble with a process zone fracture propagation criteria. The effect of the yield strength on the behavior of this curve is shown with two different examples; (a) ry ¼ 1:5  104 Nm2 and (b) ry ¼ 5  103 Nm2 . The higher yield strength (a) produces a curve where the drops in pressure during fracture events become less obvious as the bubble grows and the curve appears as a slightly decreasing flat-line pressure record. The lower yield strength produces a ‘‘saw-toothed” pressure record which is the archetypical LEFM fracture pressure record.

2590

C.K. Algar, B.P. Boudreau / Geochimica et Cosmochimica Acta 73 (2009) 2581–2591

Fig. 11. Bubble growth curves with a process-zone fracture criteria. Three different values of yield strength were examined, ry ¼ 5:0  103 Nm2 ; ry ¼ 1:0  104 Nm2 ; ry ¼ 1:5  104 Nm2 and compared to growth with a constant post-fracture pressure criteria. This shows that weaker sediments (lower yield stress) grow faster then stronger (high yield stress) sediments due to the increased length of fracture propagation.

range of parameter values that sediments will encompass, it should not be surprising that both types of behavior are common (even in the same core sample). Fig. 11 explores the sensitivity of the growth rate to the yield strength, and this is compared to growth under the constant post-fracture pressure propagation rule. From this we can see that higher yield strength values slow down growth by decreasing the distance that fractures propagate. Lower values give growth times similar to the constant post-fracture pressure criteria. The process zone propagation rule does not drastically change the growth rates of bubbles, but nevertheless over the range of yield strength values examined, the growth times of a 0.5 cm3 reference bubble varied by approximately a day, a difference of about 20%. As a result yield stress measurements on sediments should be made in the future. 8. CONCLUSIONS We have tested the assumption of quasi-steady state employed in the previous LEFM reaction–diffusion model of bubble growth and have found that the transient nature of the methane concentration field surrounding a growing bubble cannot be ignored. Furthermore, we have found that the use of the fully transient form of the reaction–diffusion model, combined with the proper form of the COD equation (Eq. 4), results in about a 40% decrease in the time required to grow a bubble to a reference size of 0.5 cm3 for CLB parameters. Bubbles of the size released at CLB will grow on daily time scales. Though growth rates of natural bubbles are scarce, agreement with the observations in Martens and Klump (1980) and Klein (2006) suggests that

our predictions are inline with what could be expected in nature. In general, this paper highlights the need for better characterization of the sediment environments in which these bubbles occur. In addition to measuring the relevant mechanical properties mentioned in this paper, attempts should eventually be made to link these properties to more commonly measured properties such as grain size or organic matter content. To date the only reported values for the mechanical properties relevant to bubble growth are reported in Johnson et al. (2002) for a single site at Cow Bay, (Cole Harbour), Nova Scotia, and we have applied these values to the CLB site. It is obvious, though, that not all sediments are the same and that there will be variability in these parameters. For this reason we performed simulations for a variety of parameter values, highlighting how different sediments could potentially influence bubble growth. Finally we have proposed a criteria for fracture propagation which is consistent with the principles of LEFM and found that the adjusted theory reproduces the range of pressure records seen in bubble growth experiments. (Johnson et al., 2002). In closing, much of our predictions are at present speculative; as time goes on and the mechanical properties of sediments become better known, the potential for more qualitative predictions of bubble growth will arise. Our model should then go a long way in helping to better understand the influence that bubbles have on their surroundings, both on a local scale (solute exchange processes) and on larger scales, such as methane fluxes to the atmosphere. ACKNOWLEDGMENTS This research was funded by the US Office of Naval research through Grants N00014-08-1-0818 and N00014-05-1-0175 (Project Managers, J. Eckman and T. Drake). Support was also provided by the Natural Sciences and Engineering Research Council of Canada and the Killam Trust. We would also like to thank Dr. C. Ruppel and Dr. F. Burdige for their insightful comments and suggestions.

REFERENCES Abegg F. and Anderson A. L. (1997) The acoustic turbid layer in muddy sediments of Ecernf} orde Bay, Western Baltic: methane concentration, saturation and bubble characteristics. Mar. Geol. 137, 137–147. Albert D. B., Martens C. S. and Alperin M. J. (1998) Biogeochemical process controlling methane in gassy coastal sediments—Part 2: groundwater flow control of acoustic turbidity in Eckernf} orde Bay sediments. Cont. Shelf Res. 18, 1771–1793. Anderson A. L. and Hampton L. D. (1980a) Acoustics of gas bearing sediments. I. Background. J. Acoust. Soc. Am. 67, 1865–1889. Anderson A. L. and Hampton L. D. (1980b) Acoustics of gas bearing sediments. II. Measurements and models. J. Acoust. Soc. Am. 67, 1890–1903. Anderson A. L., Abegg F., Hawkins J. A., Duncan M. E. and Lyons A. P. (1998) Bubble populations and acoustic interaction with the gassy floor of Eckernf} orde Bay. Cont. Shelf Res. 18, 1807–1838.

Bubble Growth in Sediments Boudreau B. P., Gardiner B. S. and Johnson B. D. (2001) Rate of growth of isolated bubbles in sediments with a diagenetic source of methane. Limnol. Oceanogr. 43(3), 616–622. Boudreau B. P., Algar C., Johnson B. D., Croudace I., Reed A., Furukawa Y., Dorgan K. M., Jumars P. A., Grader A. S. and Gardiner B. S. (2005) Bubble growth and rise in soft sediments. Geology 33(6), 517–520. Broek D. (1982) Elementary Engineering Fracture Mechanics. Narinus Nijhoff, Boston. Brook E. J., Showers T. and Orchardo J. (1996) Rapid variations in atmosphere methane concentration during the past 110,000 years. Science 273, 1087–1091. Chanton J. P., Martens C. S. and Kelley C. A. (1989) Gas transport from methane-saturated, tidal freshwater and wetland sediments. Limnol. Oceanogr. 34(5), 807–819. Crill P. M. and Martens C. S. (1986) Methane production from bicarbonate and acetate in an anoxic marine sediment. Geochim. Cosmochim. Acta 50(9), 2089–2097. Gardiner B. S., Boudreau B. P. and Johnson B. D. (2003a) Growth of disk-shaped bubbles in sediments. Geochim. Cosmochim. Acta 67(8), 1485–1494. Gardiner B. S., Boudreau B. P. and Johnson B. D. (2003b) Slow growth of an isolated disk-shaped bubble of constant eccentricity in the presence of a distributed gas source. Appl. Math. Model. 27, 817–829. Gross D. and Seelig T. (2006) Fracture Mechanics. With an Introduction to Micromechanics, Springer, Berlin. Haeckel M. P., Boudreau B. and Wallmann K. (2007) Bubbleinduced porewater mixing: a 3-D model for deep porewater irrigation. Geochim. Cosmochim. Acta 71(21), 5135–5154. Heimpel H. and Olson P. (1994) Buoyancy-driven fracture and magma transport through the lithosphere: models and experiments. In Magmatic Systems (ed. M. P. Ryan). Academic Press, San-Diego, USA, pp. 223–240. Hsieh D. Y. and Plesset M. S. (1961) Theory of rectified diffusion of mass into gas bubbles. J. Acoust. Soc. Am. 33, 206–215. Jackson J. R., Williams K. L., Wever T. F., Friedrichs C. T. and Wright L. D. (1998) Sonar evidence for methane ebullition in Eckernf} orde Bay. Cont. Shelf Res. 18, 1893–1915.

2591

Jain A. K. and Juanes R. (2008) Pore-scale mechanistic study of the preferential mode of hydrate formation in sediments: coupling of multiphase fluid flow and sediment mechanics. In Proceedings of the Sixth International Conference on Gas Hydrates (ICGH 2008), Vancouver, 6–10 July 2008. Johnson B. D., Boudreau B. P., Gardiner B. S. and Mass R. (2002) Mechanical response of sediments to bubble growth. Mar. Geol. 187, 247–363. Judd A. G. (2004) Natural seabed gas seeps as sources of atmospheric methane. Environ. Geol. 46, 988–996. Klein Stephe (2006) Sediment porewater exchange and solute release during ebullition. Mar. Chem. 102, 60–71. Lawn B. R. and Wilshaw T. R. (1975) Fracture of Brittle Solids. Cambridge University Press, Cambridge, England. Lebedev N. N. (1972) Special Functions and their Applications (ed. and trans. R. A. Silverman). Dovor Publishing, New York. Martens C. S. and Klump J. V. (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. and Goldhaber M. B. (1978) Early diagenesis in transitional sedimentary environments of the White Oak River Estuary, North Carolina. Limnol. Oceanogr. 23(3), 428–441. Rubin A. M. (1993) Tensile fracture of rock at high confining pressure: Implications for dike propagation. J. Geophys. Res. 98, 15919–15935. Sills G. C. and Wheeler S. J. (1992) The significance of gas from offshore operations. Cont. Shelf Res. 12, 1239–1250. 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. Timoshenko S. P. and Goodier J. N. (1970) Theory of Elasticity. McGraw-Hill, New York. Wever T. H., Abegg F., Fiedler H. M., Fechner G. and Stender I. H. (1998) Shallow gas in the muddy sediments of Eckernf} orde Bay, Germany. Cont. Shelf Res. 18, 1715–1739. Associate editor: David J. Burdige