Non-Fickian Diffusion and the Accumulation of Methane Bubbles in Deep-Water Sediments
In the absence of fractures, methane bubbles in deep-water sediments can be immovably trapped within a porous matrix by surface tension. The dominant mechanism of transfer of gas mass therefore becomes the diffusion of gas molecules through porewater. The accurate description of this process requires non-Fickian diffusion to be accounted for, including both thermodiffusion and gravitational action. We evaluate the diffusive flux of aqueous methane considering non-Fickian diffusion and predict the existence of extensive bubble mass accumulation zones within deep-water sediments. The limitation on the hydrate deposit capacity is revealed; too weak deposits cannot reach the base of the hydrate stability zone and form any bubbly horizon.
pacs:66.10.C-Diffusion and thermal diffusion and 47.56.+rFluid flow through porous media and 91.50.HcMarine geology: gas and hydrate systems
The occurrence of methane bubbles within porous water-saturated sediments is widespread around the ocean margins. The gas within them plays an important role in both submarine hazards, such as submarine landslides Kayen-Lee-1991 () as well as the formation of resources, such as methane-hydrate deposits Davie-Buffett-2001 (); Davie-Buffett-2003b (); Archer-2007 (). The stability of the bubbles has a significant control on the methane flux from the sediments into the ocean-atmosphere system.
In porous sediments the bubbles are trapped within the matrix pores. Large moving bubbles are unstable, as they split into smaller bubbles during migration Lyubimov-etal-2009 (); Barry-etal-2010 (), and smaller bubbles are trapped by pore-throats or by surface tension forces. The minimum pore throat diameter required to trap a small bubble, when there is no strong pumping of the fluid through the porous matrix, can be calculated. The surface tension forces ( is the surface tension) should overwhelm the buoyancy force ( is the density, is the gravity), i.e., . For gas-water systems, one finds, . Making allowance for the inhomogeneity of pores and the geometry of contacts, one should decrease this estimate to , which still suggests trapping even for sands. For a soft mud, mechanics of bubbles and porous matrix can be different Algar-Boudreau-Barry-2011 () and is not considered here.
As bubbles are trapped, and in the absence of significant groundwater movement transporting dissolved methane, the dominant mechanism of methane mass transfer in deep-water sediment is by diffusion of methane through porewater, controlled by (i) the methane saturation of the aqueous solution throughout the sediment volume, and (ii) non-Fickian diffusion laws. The latter pertain to two processes: firstly, the geothermal gradient causes thermodiffusion (the Soret effect Soret-1879 ()), where the temperature gradient induces solute flux (as recognized in several fields, e.g. Richter-1972 ()); and secondly, the impact of gravity on dissolved molecules. These non-Fickian contributions to the diffusive flux mean that the solute flux cannot solely be determined by the gradient of the solute concentration.
Currently non-Fickian diffusion is rarely considered in the modelling of deep-water sediment systems (Davie-Buffett-2001 (); Davie-Buffett-2003b (); Haacke-Westbrook-Riley-2008 (); Garg_etal-2008 (), and others involved in gas hydrate modelling only address Fickian diffusion). However, we suggest that non-Fickian processes are important; in particular, they may cause methane to migrate against the direction of the steepest decrease of concentration under certain conditions. Importance of non-Fickian diffusion was also demonstrated for the salinity transport in hydrate-bearing sediments Goldobin-CRM-2013 (). Unfortunately, the experimental value of the thermodiffusion coefficient for the aqueous methane solution is unknown and can only be roughly assessed theoretically. Therefore, we treat it as a free parameter in our investigation. This uncertainty also justifies the fact that prior numerical modelling work did not include non-Fickian processes.
In researching the horizontal through-flow of water and vertical aqueous oxygen transport in porous bubble-bearing sediments, Donaldson-etal-1997 () considered the hydrodynamic dispersion (or “turbulent diffusion” Bernard-Wilhelm-1950 ()), caused by water transport through irregular pore channels, as a transport mechanism and reasonably neglected the molecular diffusion. Indeed, the turbulent diffusion plays a significant role in vertical dissolved gas transport near the earth’s surface. However, it becomes insignificant away from the sediment-water interface in deep-water sediments, where vertical and horizontal displacements of water are comparable. Here it is significantly smaller than molecular diffusion—even in sandy sediments.
In the present study we consider the diffusive migration of methane in sea-floor sediments where water is saturated with methane, and some methane is gaseous (forming bubbles; see the sketch in Fig. 1a). In deep-water sediments the bubbly zone is overlaid by the methane hydrate stability zone, where non-dissolved methane forms not gaseous bubbles but hydrate, which is included in our treatment as well.
2 Transport processes in sediments
On the field scale, deep-water sediments are typically much more uniform horizontally rather than vertically. Consequently, we consider a system that is uniform horizontally. The depth below the water-sediment interface is measured by the -coordinate (Fig. 1).
2.1 Diffusion in non-isothermal aqueous solutions
Here is the mass fraction of the solute in the solvent, is the solute molecular diffusion coefficient, is the porosity of the solid matrix, is the tortuosity factor featuring the pore geometry, is the fluid density. The first term describes the “ordinary” Fickian diffusion, . The second term represents the thermodiffusion effect appearing in non-isothermal systems, where temperature inhomogeneity causes a solute flux. The strength of the thermodiffusion effect is characterized by the thermodiffusion constant (the conventional Soret or separation coefficient ). The third term describes the action of gravity on solute molecules; is the universal gas constant, , and are the molar masses of the solute and solvent, respectively, and is the number of solvent molecules in the volume occupied by one solute molecule in the solution. The value of can be precisely derived for – systems from the dependence of the solution density on its concentration Hnedkovsky-Wood-Majer-1996 (); one obtains and .
When the liquid is saturated with gas bubbles, the concentration of solute in solvent equals the solubility, , throughout the liquid volume in the bubble-bearing zone because the bubbles are in local thermodynamic equilibrium with the solution. Thus, the solute flux depends merely on the temperature and pressure fields, and , and the solution concentration is not a free variable, . At high pressure and low enough temperature, the hydrate form is more thermodynamically preferable for methane than the gaseous form. The hydrate zone can overlie the bubbly-bearing zone (Fig. 1). In the presence of hydrate the solution concentration equals the solubility at equilibrium with hydrate. For the calculation of the solubility (in bubbly and hydrate zones), half-empiric models developed in Duan-Mao-2006 (); Sun-Duan-2007 () are employed.
2.2 Pressure and temperature
The solubility profiles and thus transport processes depend on the pressure and temperature profiles. The system is essentially characterized by the hydrostatic pressure and geothermal temperature gradient . Although the role of porosity nonuniformity for the geothermal gradient was demonstrated in Goldobin-EPL-2011 (), we follow the conventional approximation of a linear temperature profile Davie-Buffett-2001 (); Garg_etal-2008 (); Haacke-Westbrook-Riley-2008 (). Therefore,
where is the atmospheric pressure, is the depth of the water body above the bubble-bearing porous sediments, and is the temperature at the sediment-water interface.
2.3 Sediment compaction and transport processes
The sedimentation process and the weight of the above-laying sediments result in a non-uniform porosity profile and a non-uniform velocity of the downward sediment motion. The sediment porosity is typically adopted in the form Davie-Buffett-2001 (); Haacke-Westbrook-Riley-2008 (), and the resulting sediment velocity is Davie-Buffett-2001 ()
where is the sedimentation rate.
The sedimentation and compaction processes create an ascending filtration flux of water through the sediments. The net water mass flux is contributed by the filtration flux and by the water transport in hydrate, which moves with sediments;
where is the hydrate density, is the mass fraction of in methane hydrate, and is the volumetric fraction of hydrate in pores or hydrate saturation.
Whilst the volumetric fraction of bubbles (or gas saturation) in pores, say , is vanishingly small, bubbles are immovably trapped in pores and move with sediments Abaci-Edwards-Whittaker-1992 (); Davie-Buffett-2001 (); Haacke-Westbrook-Riley-2008 (). However, transport processes in the system can result in a gas saturation which exceeds the critical value . When the critical saturation is exceeded, gas leakage occurs.
The net methane mass flux is contributed by the ascending filtration flux of aqueous solution, the motion of gas bubbles and/or hydrate with sediments, and molecular diffusion in aqueous solution. In the absence of gas leakage, it reads
where the tortuosity factor is assumed Haacke-Westbrook-Riley-2008 (), is the mass fraction of in hydrate, and
Eq. (5) with corresponds to the Fickian diffusion law, , and characterizes the strength of non-Fickian part of the solute flux. As the value of is unknown, is treated as a free parameter in our study. In our model, the zones of hydrate and gas bubbles do not overlap; either or can be non-zero at a given location.
While the diffusion coefficient is nearly independent of pressure Sachs-1998 (), its dependence on temperature is not to be neglected as the temperature change from to causes the diffusivity increase by factor . The formula with , , , well fits experimental data summarized by Sachs-1998 ().
In the present work we do not consider details of the process of methane generation from the organic part of sediments, assuming that it takes place in the upper part of sediments and methane is simply present at certain depth of about below the water-sediment interface (cf Haacke-Westbrook-Riley-2008 ()).
Both the numerical modelling with finite difference method and analytical treatment were performed. The results of numerical simulation (also including different gas leakage models) are in agreement with the analytical theory (as can be later seen from Fig. 2d). It is convenient to start the presentation of results with time-independent solutions to the problem.
|-folding depth of porosity|
|fluid filtration velocity|
3.1 No-leakage time-independent states
Three zones with different transport features can be distinguished
(see Fig. 1):
(1) the lowermost zone of undersaturated aqueous solution of methane without gas bubbles, ,
(2) the zone of saturated solution bearing gas bubbles, , and
(3) the lower part of hydrate stability zone (HSZ) with saturated solution and hydrate, .
In the lower part of HSZ we consider, no process of conversion of organic sediments into methane occurs; we assume the conversion to be finished above this zone. In Fig. 1b, one can see a sample profile of the net mass fraction of methane in pores composed of dissolved methane and methane in hydrate or gas bubbles.
At the time-independent state, neither the hydrate distribution nor the free-gas one change with time, which requires the methane mass flux (5) to be uniform along the sediment column. As we do not consider methane influx from deep massifs, i.e., methane mass flux is zero deep in sediments, this uniform flux should be zero everywhere. The water flux (4) is uniform, , where is the filtration velocity below HSZ (cf Davie-Buffett-2001 ()).
For zero methane flux and given water flux, one can find the solution concentration in zone (1);
In zone (2), the solute concentration equals the solubility, , and, in the absence of leakage, gas saturation
In zone (3), the solute concentration equals the solubility, , and hydrate saturation
(Eqs. (8) and (9) are derived from Eqs. (4) and (5) for small and .) In Fig. 1b the time-independent state is plotted for , the graphs for different are qualitatively similar. Gas saturation immediately above the lower boundary of the bubbly zone is finite and enough for the formation of the second bottom simulating reflector of seismic waves, which is not associated with hydrate.
3.2 Evolution of the bubbly zone and persistence of marine hydrate deposits
For a time-independent state determined by Eqs. (7)–(9) the position of the boundary between the bubbly zone and the zone of undersaturated solution is not imposed; all three profile plotted is Fig. 2a satisfy the equations. Now let us consider small deviation from these profiles.
When hydrate saturation at the base of HSZ is smaller than that in the time-independent state, upward methane flux appears. Indeed, the diffusion flux is the same as in the time-independent state, because it is controlled by the solubility profile, while we have strictly downward transfer of hydrate with sediments, and this transfer is diminished in comparison to the transfer in the time-independent state. In the time-independent state these two fluxes are balanced to zero net flux, while with the deficiency of hydrate we have the deficiency of downward flux and thus the net flux is upward. An upward net flux of methane results in depletion of methane deep in sediments and gradual retreat of the bubbly zone. After some period of time the bubbly zone disappears and, moreover, hydrate zone retreats from the base of HSZ to the area of methane generation from sediments. In this case one observes HSZ with undersaturated aqueous solution and no hydrate at the base of it and no bubbly zone; this has been reported to be widespread in nature Proc_ODP-V164 (). Thus, one can see that the hydrate deposit cannot be too weak, otherwise it is not persistent. In Fig. 3 the minimal hydrate saturation is plotted as a function of the water-body depth; this value is also affected by the non-Fickian drift strength , for a larger a stronger deposit is required.
When the hydrate saturation at the base of HSZ exceeds the critical hydrate saturation, one observes an opposite situation: a downward net flux of methane and a bubbly zone gradually advancing into deep sediments. However, the advance of the bubbly zone is not unbounded. In Fig. 2b, one can see that for time-independent states the gas saturation increases with depth and at certain depth, , exceeds the critical gas saturation. The critical gas saturation, above which gas leakage starts, depends on many factors, including the pore geometry and the solid matrix material Abaci-Edwards-Whittaker-1992 (), and the rate of gas release from the solution (for discussion see Haacke-Westbrook-Riley-2008 ()). The critical gas saturation specific for different porous massifs is about (e.g., see review in Haacke-Westbrook-Riley-2008 ()), and we adopt this value for our analysis. Leakage is a hydrodynamical transport and as such is much more efficient than the molecular diffusion transport or the motion with sediments. Therefore, any significant excess of gas saturation over the critical value results in a mass flux which cannot be balances by molecular diffusion or sedimentation, and, on the time scales of the hydrate deposit formation, gas saturation can be only under- or nearly-critical. The bubbly zone cannot advance deeper, as a larger gas saturation is required to have a downward methane transfer below depth . Leaking gas is accumulated just above the leakage zone until the critical saturation is reached there. Then leakage zone extends, and finally entire bubbly zone becomes the leakage zone with gas saturation slightly exceeding the critical saturation, as shown in Fig. 2c. Hence, for persistent hydrate deposits () one can observe the formation of extensive bubbly zone with critical gas saturation; the lower boundary of this zone is determined by condition
In Fig. 2d, one can compare the analytical time-independent profile to the results of numerical simulation with a comprehensive account of processes: the rate of methane generation by anaerobic bacteria
with and (cf Davie-Buffett-2001 ()) and the gas leakage within the bubbly zone with relative permeability for the gas phase Mavko-Mukerji-Dvorkin-2009 (). The numerically calculated profile slightly deviates from the analytical one for realistic sediment permeability and this deviation remains non-large even for as small permeability as . Noteworthy, the base of the free-gas zone remains unshifted even for small permeability.
In numerical simulation the described evolution was observed with different leakage laws, which may strongly vary from system to system (e.g. Abaci-Edwards-Whittaker-1992 ()), and realistic values of massif permeability. As explained and demonstrated in Fig. 2d, the bubbly zone is tolerant to specific features of the leakage law with given ; realistic model parameters can only change the small excess of the gas saturation over the critical value. In our analytical theory we do not discuss the narrow hydrate–gas recycling zone immediately below HSZ, where gas saturation may significantly exceed the critical saturation. Features of this zone depend on the leakage model which is highly uncertain without thorough knowledge of the massif properties and is site-specific.
In Fig. 4, one can see that the location of the extensive free-gas (bubbly) zone is affected by the non-Fickian drift strength . Moreover, the minimal depth of the water body above sediments needed for the extensive bubbly zone to appear is also controlled by . An extensive bubbly zone is not possible below HSZ beneath shallow water bodies because the time-independent state requires gas saturation higher than the critical saturation and cannot be maintained because of gas leakage. It is noteworthy that, due to hydrate–gas recycling, a narrow bubbly layer immediately below HSZ will be still presented wherever hydrate is present at the base of HSZ.
3.3 Importance and uncertainties
The behavior described is considerably influenced by the non-Fickian drift of methane. We suggest that the Fickian diffusion law, which has been adopted in Davie-Buffett-2001 () and its successors and corresponds to , should be modified.
The unresolved issue here is the specific value of for methane. The authors are not aware of experimental data on thermodiffusion of methane in water, though there are a lot of experimental studies on the thermodiffusion of methane in mixtures of hydrocarbons (e.g. Wittko-Koehler-2005 ()). Theoretical studies (e.g. Semenov-2010 ()) provide formulae for calculation of the thermodiffusion constant from inter-molecular potentials which are poorly established for water because of hydrogen bonds. We can only calculate (for ) and use a rough conjecture that one can expect Goldobin-Brilliantov-2011 (). Hence, can be expected for aqueous solutions with geothermal gradient .
We have theoretically explored the process of diffusive migration of aqueous methane in the presence of bubbles, when they are immovably trapped by a porous matrix—as occurs commonly in seafloor sediments, swamps, or terrestrial aquifers. The effect of temperature inhomogeneity across the system (geothermal gradient) and gravitational force have been accounted for.
Non-Fickian corrections—thermodiffusion and gravitational segregation—appear to play an important role in the migration of methane in sediments in deep-water settings. The positive thermodiffusion effect (, cf Eq. (6)) and the gravitations segregation of methane in water make negative contribution to , and, therefore, they assist the formation of an extensive methane gas accumulation zone in the upper part of the sediment column under deep water bodies. For instance, Fig. 4 illustrates that, for conditions of the west Svalbard continental slope Haacke-Westbrook-Riley-2008 (), non-Fickian diffusion can either extend () or shrink () the zone of methane gas accumulation.
Remarkably, hydrate deposits with too small hydrate saturation at the base of the hydrate stability zone should suffer diffusive depletion and retreat from the base of HSZ to the region of methane generation from sediments or completely disappear. This explains why some natural hydrate deposits are reported to possess no hydrate and no bottom simulating reflector at the base of HSZ. The positive thermodiffusion effect and gravitational segregation, both resulting in negative , decrease the minimally required hydrate saturation, as can be seen in Fig. 3.
Unfortunately, we cannot determine precise values for the thermodiffusion constant of aqueous solutions of methane from the literature, and can only rely on theoretical predictions (e.g. Semenov-2010 (); Goldobin-Brilliantov-2011 ()), to estimate their values, as we do here. Our findings highlight the necessity of experimental determinations of the thermodiffusion constant for aqueous methane solutions.
Acknowledgements.The work has been supported by NERC (NE/F021941/1). DSG acknowledges the financial support by the Government of Perm Region (Contract C-26/212). This paper is published with the permission of the Director of the British Geological Survey.
- (1) R. E. Kayen and H. J. Lee, Pleistocene slope instability of gas hydrate-laden sediment on the Beaufort Sea margin, Marine Geotechnology 10, (1991) 125–141.
- (2) M. K. Davie and B. A. Buffett, A numerical model for the formation of gas hydrate below the seafloor, J. Geophys. Res. 106, (2001) 497–514.
- (3) M. K. Davie and B. A. Buffett, Sources of methane for marine gas hydrate: inferences from a comparison of observations and numerical models, Earth Planet. Sci. Lett. 206, (2003) 51–63.
- (4) D. Archer, Methane hydrate stability and anthropogenic climate change, Biogeosciences 4, (2007) 521–544.
- (5) D. V. Lyubimov, S. Shklyaev, T. P. Lyubimova and O. Zikanov, Instability of a drop moving in a Brinkman porous medium, Phys. Fluids 21, (2009) 014105.
- (6) M. A. Barry, B. P. Boudreau, B. D. Johnson and A. H. Reed, First-order description of the mechanical fracture behavior of fine-grained surficial marine sediments during gas bubble growth, J. Geophys. Res. 115, (2010) F04029.
- (7) C. K. Algar, B. P. Boudreau and M. A. Barry, Initial rise of bubbles in cohesive sediments by a process of viscoelastic fracture, J. Geophys. Res. 116 (2011) B04207.
- (8) C. Sorét, Sur l’état d’équilibre que prend, au point de vue de sa concentration, une dissolution saaline primitivement homogène, dont deux parties sont portées à des températures différentes, Archives des Sciences Physiques et Naturelles de Genève 2, (1879) 48–61.
- (9) J. Richter, Evidence for significance of other-than-normal diffusion transport in soil gas exchange, Geoderma 8, (1972) 95–101.
- (10) R. R. Haacke, G. K. Westbrook and M. S. Riley, Controls on the formation and stability of gas hydrate-related bottom-simulating reflectors (BSRs): A case study from the west Svalbard continental slope, J. Geophys. Res. 113, (2008) B05104.
- (11) S. K. Garg, J. W. Pritchett, A. Katoh, K. Baba and T. Fujii, A mathematical model for the formation and dissociation of methane hydrates in the marine environment, J. Geophys. Res. 113, (2008) B01201.
- (12) C. Sorét, Sur l’état d’équilibre que prend, au point de vue de sa concentration, une dissolution saaline primitivement homogène, dont deux parties sont portées à des températures différentes, Archives des Sciences Physiques et Naturelles de Genève 2, (1879) 48–61.
- (13) J. Richter, Evidence for significance of other-than-normal diffusion transport in soil gas exchange, Geoderma 8, (1972) 95–101.
- (14) R. R. Haacke, G. K. Westbrook and M. S. Riley, Controls on the formation and stability of gas hydrate-related bottom-simulating reflectors (BSRs): A case study from the west Svalbard continental slope, J. Geophys. Res. 113, (2008) B05104.
- (15) S. K. Garg, J. W. Pritchett, A. Katoh, K. Baba and T. Fujii, A mathematical model for the formation and dissociation of methane hydrates in the marine environment, J. Geophys. Res. 113, (2008) B01201.
- (16) D. S. Goldobin, Non-Fickian diffusion affects the relation between the salinity and hydrate capacity profiles in marine sediments, Compte Rendus Mecanique 341, (2013) 386–392.
- (17) J. H. Donaldson, J. D. Istok, M. D. Humphrey, K. T. O’Reilly, C. A. Hawelka and D. H. Mohr, Development and testing of a kinetic model for oxygen transport in porous media in the presence of trapped gas, Ground Water 35, (1997) 270–279.
- (18) R. A. Bernard and R. H. Wilhelm, Turbulent diffusion in fixed beds of packed solids, Chem. Eng. Prog. 46, (1950) 233–244.
- (19) R. B. Bird, , W. E. Stewart and E. N. Lightfoot, Transport Phenomena (Wiley, 2007).
- (20) L. Hnedkovsky, , R. H. Wood and V. Majer, Volumes of aqueous solutions of CH4, CO2, H2S and NH3 at temperatures from 298.15 K to 705 K and pressures to 35 MPa, J. Chem. Thermodyn. 28, (1996) 125–142.
- (21) Z. Duan and S. Mao, A thermodynamic model for calculating methane solubility, density and gas phase composition of methane-bearing aqueous fluids from 273 to 523 K and from 1 to 2000 bar, Geochim. Cosmochim. Acta 70, (2006) 3369–3386.
- (22) R. Sun and Z. Duan, An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments, Chem. Geol. 244, (2007) 248–262.
- (23) D. S. Goldobin, Scaling of transport coefficients of porous media under compaction, Europhys. Lett. 95, (2011) 64004.
- (24) S. Abaci, J. Edwards and B. Whittaker, Relative permeability measurements for two phase flow in unconsolidated sands, Mine Water Environ. 11, (1992) 11–26.
- (25) W. Sachs, The diffusional transport of methane in liquid water: method and result of experimental investigation at elevated pressure, J. Pet. Sci. Eng. 21, (1998) 153–164.
- (26) C. K. Paull, R. Matsumoto, P. J. Wallace and W. P. Dillon (Eds.), Proceedings of the Ocean Drilling Program, in: Scientific Results, vol. 164 (Ocean Drilling Program, College Station, TX, 2000).
- (27) G. Mavko, T. Mukerji and J. Dvorkin, The Rock Physics Handbook (Cambridge University Press, 2009).
- (28) G. Wittko and W. Köhler, Universal isotope effect in thermal diffusion of mixtures containing cyclohexane and cyclohexane-d[sub 12], J. Chem. Phys. 123, (2005) 014506.
- (29) S. N. Semenov, Statistical thermodynamic expression for the Soret coefficient, Europhys. Lett. 90, (2010) 56002.
- (30) D. S. Goldobin and N. V. Brilliantov, Diffusive migration of mass in bubbly media, Phys. Rev. E 84, (2011) 056328.