On the fate of a CO

On the fate of a CO

lake disposed in the deep ocean

Technical report prepared by Dr. Ilker Fer

as a contribution to the COSMOS (Study on the development of new CO2 sending method for Ocean Storage), subproject on hydrodynamic

and oceanographic modeling

Geophysical Institute, University of Bergen March 2002

1. Introduction

Enhanced emission of greenhouse gases, particularly carbon dioxide (CO2), to the atmosphere is widely accepted to affect the global climate system [Houghton et al., 1995]. The atmospheric CO2 content at present is about 25% higher than pre-industrial levels. Over the past two decades, multidisciplinary research has been intensified to stabilize the CO2 level in the atmosphere. One of the potential options to mitigate atmospheric CO2 level is to capture CO2

from fossil fuel combustors and to purposefully dispose and sequester elsewhere (e.g., in ocean, deep saline aquifers, depleted gas and oil wells, and coal beds, etc.). Ocean appears to be the most preferable option being the largest potential sink for anthropogenic CO2 and is,

particularly the deep ocean, highly unsaturated in CO2. Marchetti [1977] was the first to propose the ocean disposal of CO2 in order to accelerate the natural ocean uptake of

atmospheric CO2. He suggested that an efficient long-term sequestration could be achieved through the Gibraltar Strait where the outflow of dense water cascades to ∼1000 m depth which, in consequence, spreads out in the North Atlantic.

The research on ocean disposal option have mostly focused on predicting the behavior and the dissolution time scale of the released CO2 [see e.g., Handa and Ohsumi, 1995] and on quantifying the environmental impacts to marine systems [Adams et al., 1997; Tamburri et al., 2000]. Different scenarios of CO2 disposal in ocean have been proposed at various depths and in different forms in relation to the phase properties of CO2. The phase diagram for the CO2 – water system shows that when pressure is greater than 44.5 bar and the temperature is less than 9.85 °C, clathrate-hydrate crystal (hydrate, hereinafter) develops. Throughout this report, pressure is given in units of bar, atm, or MPa (1 bar = 0.9869 atm = 0.1 MPa) whichever is convenient, particularly with respect to reported results by other researchers. Density profiles of liquid CO2, CO2-saturated seawater, and seawater are shown in Figure 1, together with

approximate density of hydrate, condensation depth (uppermost arrow on the pressure axis), and the depth at which we might expect hydrate to develop (lowermost arrow on the pressure axis). Because of its large compressibility, liquid CO2 is relatively denser than seawater at depths greater than ∼3000 m. Liquefied CO2 released at shallower depths will lead to rising droplet plumes which upon reaching the condensation depth of ∼ 450 m will form CO2 bubbles [Liro et al., 1992]. Combined effect of rising droplets and sinking, dense, CO2–enriched

seawater has been recently simulated by Alendal and Drange [2001]. Furthermore, dissolution of CO2 increases the density of water and it was suggested by Drange et al. [1993] that a gravity current cold be achieved towards the deep ocean if sufficiently dense CO2 enriched water [Haugan and Drange, 1992] were released on a sloping bottom. Disposal of liquid CO2 at depths where it is denser than seawater (> 3000 m) is expected to fill topographic depressions, which in turn accumulates as a large lake of CO2 [Ohsumi, 1993; Nakashiki, 1997] over which a thin hydrate layer forms and retards the dissolution. The hydrate films are not perfect insulators against the inter-liquid-phase CO2 transfer and they are continuously ‘metabolized’ through decomposition of aged crystals and formation of new crystals [Mori, 1998]. The dissolved CO2

diffuses out into the benthic boundary layer, where advection and mixing takes place due to bottom current and turbulence near the bottom.

liquid CO₂ 0 ^oC 10 ^oC

850 900 950 1000 1050 1100 1150

seawater CO₂-saturated seawater 0

100

200

300

400

0

1000

2000

3000

4000

Depth (m)

Pressure (bar)

Density (kg m^-3)

hydrate

Figure 1. Density profiles of liquid CO2 at 0°C and 10°C (solid curves), CO2–saturated seawater (dashed line), and seawater (dotted line). Densities of seawater and CO2-saturated seawater are

calculated for T = 5°C and S = 35 psu using the international equation of state [UNESCO, 1981], and the relation (25) which incorporates the effect of dissolved CO₂ concentration [Weiss, 1974; Giggenbach, 1990], respectively. Those for the liquid CO₂ are given by Vukalovich and Altunin [1968]. Approximate density of hydrate [Udachin et al., 2001] is shown by the dot on the density axis. The phase boundary for gas-liquid CO₂ transition, i.e. condensation pressure, [Lide, 2000] and for hydrate formation [King, 1969] at 5°C are shown by the uppermost and lowermost arrows on the pressure axis, respectively.

Depth is calculated from pressure.

In this report, we address the option of CO2 disposal as a lake in the deep ocean. We attempt to envisage the fate of the liquid CO2 pool through simple numerical simulations of the two-dimensional advection-diffusion equation, incorporating the bottom boundary layer

dynamics and the effects of the hydrate layer. In section 2, relevant literature is reviewed with respect to CO2 – seawater system and the boundary layer turbulence. The model is described in section 3. Subsequently, results are summarized in section 4 preceding the discussion given in section 5.

2. Background

2.1. CO2 – sea water system

Hydrate formation at liquid CO₂ – water interface

CO2 hydrate, a binary clathrate compound, forms when the cavities of hydrogen-bonded structures are occupied by CO2 molecules. The chemical formula yields CO2·5.75·H2O if all the cavities are occupied, where n = 5.75 is the hydration number. Several researchers used the maximum value (1130 kg m^-3) of the hydrate density, which is estimated from the perfect crystallographic lattice formation and the full occupancy of CO2 molecules into clathrate cages.

However, the CO2 hydrate formed under arbitrary conditions may yield various densities [Uchida, 1997]. Udachin et al. [2001] recently estimated the density of CO2 hydrate as 1120 ± 10 kg m^-3 for a calculated lattice constant at 4 °C and a composition of CO2·6.20·H2O. Here, the

main estimated uncertainty arises from the cage occupancies. At depths considered for direct disposal of CO2 in deep ocean, ∼3000 m, the hydrostatic pressures will be as high as ∼305 bar and the composition will tend toward complete filling of the cages [Udachin et al., 2001]. For the stability of a CO2 lake a better knowledge of the buoyancy behavior of the CO2 hydrate, hence its density, and of its mechanical properties are required.

Deep-sea disposal of CO2 is certain to lead to CO2 hydrate formation [Brewer et al., 1999]. Naturally occurring CO2 clathrate hydrates were observed by Sakai et al. [1990] during a submersible study. CO2 rich fluid bubbles emerging from the sea floor at 1335 to 1550 m depth in the mid-Okinawa Trough immediately formed hydrates upon contact with 3.8 °C seawater, which in turn developed fragile hydrate tubes standing on the sediments. The observations were followed by in situ experiments, notably by Honda et al. [1995] who observed the changes of initially solid CO2 contained in a cylinder while they descend it to 3073 m depth by a manned submersible and those by Brewer et al [1999] who reported intermittent overflows of liquid CO2 and massive hydrate formation from a beaker, initially half full with CO2, disposed to the ocean floor at 3627 m depth. Brewer et al. [1999] observed that the pool of liquid CO2 on the seafloor expanded in volume more than four times, forming hydrate, which eventually dissolved. Recently the dissolution rate of liquid CO2 which is accumulated into a small hole made by a push core on the ocean floor at 3600 m depth and 1.55 °C was directly measured by Dunk et al. [2002] using a ROV and a pH electrode. Their estimate of 1.7 µmol cm^-2 s^-1 is comparable to the dissolution rate of 3 µmol cm^-2 s^-1 determined from the shrinking rate of a rising stream of CO2 droplets at 800 m depth [Peltzer et al., 2000]. On a subsequent note, Peltzer et al. [2002] report on a massive ‘frost-heave’ which formed at one of the corrals at 3604 m, whereas a similar corral nearby showed no such behavior.

Various aspects of the structure and composition of the hydrate are studied by

spectroscopic and diffraction measurements [Uchida, 1997; Ikeda et al., 1999; Udachin et al., 2001] and the modeling of clathrate hydrate formation at the interface between liquid CO2 and water is extensively reviewed by Mori [1998]. Mori, in his review, classified hydrate film models into three different groups: diffusive suspension layer models [e.g. Shindo et al., 1995;

Teng et al., 1995], permeable solid plate models [e.g. Teng et al., 1996b; Mori and Mochizuki, 1997], and sediment particle aggregate layer model [Inoue et al., 1996]. Laboratory

experiments simulating CO2 disposal in deep ocean showed that the hydrate phase which grows at the interface between liquid CO2 and water is very thin (order of µm) but rather stable

[Shindo et al., 1995]. It was also reported that the ‘hydrate film’ never prevents but only retards the interphase mass transfer of CO2 provided that the water phase is not saturated with CO2. The dissolution mechanism being retarded through the hydrate film, so-called ‘barrier effect’, has been discussed by Mori and Mochizuki [1998]

Condensation and hydrate pressure

The condensation pressure, Pcond, as a function of temperature, T, for pure CO2 can be estimated by the second order polynomial curve fit to the data points of Lide [2000]

2 cond(T) 34.8649 0.90485T 0.0108504T

P = + + (1)

for -1.7 < T < 20 °C. Here, T and Pcond are in units of °C and bar, respectively. (1) gives the phase boundary for the gas-liquid CO2 transition. The condensation line as function of temperature and pressure is shown in Figure 2 by the dashed line.

The phase boundary for CO2 hydrate formation in freshwater can be obtained in a pressure-temperature diagram by linearly interpolating the dissociation pressures of 1.01, 12.64 and 45.09 bar at -24, 0, 10.2 °C temperatures, respectively [King, 1969]. A rough boundary for hydrate formation in fresh water is depicted in Figure 2 by the hatched area. The hydrate formation in seawater has been measured in laboratory conditions in the range of 1.6 – 10.6 °C under pressures up to 50 bar [Mitsubishi, 1990; Ohgaki et al., 1993] and are also shown in Figure 2.

Solubility of Carbon dioxide in seawater

In order to estimate the fate of purposefully discharged or stored CO2 in the ocean calculation of its solubility in seawater is necessary. In general, the solubility, Cs, of CO2 in seawater is a function of both temperature and pressure.

The solubility of CO2–gas in seawater at low pressures has been investigated

experimentally [Murray and Riley, 1971; Weiss, 1974]. In order to account for the non-ideality of CO2 it was proposed to employ a modified Henry’s law by expressing the activity of the solute in the gas phase by its fugacity, f, and by correcting the thermodynamic expansion of the solution by the dissolved gas. King [1969] presents a thorough discussion of Henry’s law and its application to real gases over a wide range of pressures. Modified Henry’s law can be used with confidence in systems up to 500 atm. According to this law the solubility of CO2 in liquids with elevated pressures can be represented as

( )

[

]

K

C_s = ₀f ₀ − ν (2)

where K0 is known as the modified Henry’s law constant or solubility coefficient, f is the fugacity, P is the pressure, P0 = 1 atm is the atmospheric pressure, ν is the partial molar volume of CO2, R is the gas constant, and T is the absolute temperature.

The solubility coefficient, K0, depends on the temperature and salinity according to the equation [Weiss, 1974]

( ) ( )

( ) ( )

[

]

3 2

1 0

100 / T B 100 / T B B S

100 / T ln A T / 100 A A K ln

+ +

+

+ +

= (3)

where K0, T, and S are in units of mol l^-1 atm^-1, °K, and psu, respectively. The A’s and B’s are constants and are summarized in Table 1.

Table 1. Constants for calculation of the solubility coefficient K₀ according to (3) in units of mol l^-1 atm^-

1 [Weiss, 1974].

A1 A2 A3 B1 B2 B3

-58.0931 90.5069 22.2940 0.027766 -0.025888 0.0050578

Theoretically, the equation of state of any gas may be written in virial form V ...

D V

C V 1 B RT PV

3

2 + +

+ +

= (4)

where B,C, and D are the second, third and fourth virial coefficients, respectively. The first term in (4) corresponds to the behavior of ideal gasses; the second term describes the action between two molecules, whereas the third one considers the interaction for three molecules and so on. At high pressures the interaction of several molecules becomes important and higher virial coefficients have to be considered.

The fugacity of a one-component non-ideal gas is given by [Levine, 1995]

( )

∫

/

ln f (5)

for constant temperature T. An expression for the fugacity which is valid for both the gas and liquid phases can be given by the exact thermodynamic relationship [Teng et al., 1996a]

( )

⎠

⎜ ⎞

⎝ + ⎛

⎟⎠

⎜ ⎞

⎝

⎛ −

⎟ +

⎠

⎜ ⎞

⎝

⎛ −

= _RT¹

∫

P /

ln ^V

f 0 (6)

Calculating f using the virial equation of state is valid to within 0.1% for pressures up to

∼10 atm. At higher pressures CO2 fugacities calculated from the standard virial equation of state become unreliable even if the calculation was carried out up to the third virial coefficient. A more sophisticated Benedict-Webb-Rubin (BWR hereinafter) equation of state [Benedict et al., 1940; 1942] allows for calculations up to at least 500 atm. The BWR equation can be written as

⎟⎠

⎜ ⎞

⎝

⎛ γ

⎟ −

⎠

⎜ ⎞

⎝

⎛ γ

+ α +

− +

− + + −

= ^{0 0}₂ ^{0 −2} _{3 6 2 3 2 2}

exp V 1 V

V T

c V

a V

a bRT V

T C A RT B V

P RT (7)

Here P (atm) is the pressure, V (l mol^-1) is the molar volume, T (°K) is the absolute

temperature, and R (l K^-1 atm mol^-1) is the gas constant. The numerical values of the constants in (7) are given in

Table 2 for pure CO2 [Bishnoi and Robinson, 1971; Weiss, 1974].

Table 2. The constants in the BWR equation of state for pure CO₂.

Constant Numerical value

A0 (l² atm mol^-2) 1.952194

B0 (l mol) 3.306517 x 10^-2

C0 (l² K² atm mol^-2) 1.704495 x 10⁵ a l³ atm mol^-3 2.526443 x 10^-1 b l² mol^-2 6.604129 x 10^-3 c l³ K² atm mol^-3 1.959201 x 10⁴ α l³ mol^-3 4.711701 x 10^-5 γ l² mol^-2 4.341442 x 10^-3 R l K^-1 atm mol^-1 8.205601 x 10^-2

Substituting from (7) for (RT/V – P) and evaluating the integral in (6) yields

( )

RT lnPV RT 1

PV exp V

2 V RT 2

c RT

c

RTV 5

a RTV

2 a bRT RTV

T C A RT P B

/ ln

2 2

3 3

5 2

2 0 0 0

−

−

⎟+

⎠

⎜ ⎞

⎝

⎛ γ

⎟ −

⎠

⎜ ⎞

⎝

⎛ γ

γ + γ −

+

+ α + −

−

= − ⁻

f

(8)

Solving (7) for the molar volume V using the known values of the system pressure, P, and the temperature, T, we can calculate the fugacity f from (8) for both the gas and the liquid phase of CO2. In calculation of f in the gas phase the system pressure, P, has to be replaced by the condensation pressure, Pcond, for pure CO2. With analogy to the modified Henry’s law, (2), the solubility of CO2 in seawater can be calculated for the gas phase as

( )

[

]

(T,P) ,S) (T,P K

C_s = _{0 0} f ₀ − (9)

and for the liquid phase as

( )

[

]

) (T,P ,S) (T,P K

C_s = _{0 0} f _{cond 0} − _cond (10)

Here, P0 = 1 atm is the atmospheric pressure. The condensation pressure, Pcond, in units of bar can be estimated by (1) and should be converted to atm for consistency.

The solubility of CO2 in seawater is calculated using (9) for the gas phase and (10) for the liquid phase and is shown in Figure 2. The calculated values are converted to mole fraction (C~_s

denotes the mole fraction of CO2 in the entirety of a given CO2 – water system) and are compared to those given by Stewart & Munjal [1970] in Figure 3. Extensive data of solubility have been obtained by many researchers over years. The data cover wide pressure and

temperature ranges, however, even at pressures favoring hydrate formation, were obtained in the absence of hydrate. Solubility data of CO2 in seawater coexisting with hydrate phase were reported only recently, notably by Aya et al. [1997]. Solubility data measured by them at 30 MPa at the presence of hydrate are also shown in Figure 3.

Fugacities calculated from (8) for 1, 10, 20 and 30 atm pressure at temperature range between 0 – 20 °C are given in Figure 4. Calculated values agree well with the data given by King [1969].

Temperature (^oC)

Pressure (bar)

200 400 600

800 1200 1000 1400

1600

2 0 2 4 6 8 10 12 14 16 18 20

1 50 100 150 200 250 300 350

2 0 2 4 6 8 10 12 14 16 18 20

1 10 20 30 40 50

400 600

800 1000

1200

200 1400

Temperature (^oC)

Pressure (bar)

Figure 2. Solubility of CO₂ in seawater with salinity of 35 psu. Solubility in gas and liquid phase are calculated using (9) and (10), respectively. Contours are in mol m^-3. The pressure scale is in bar = 10⁵ Pa

= 0.98692 atm. The dashed line is the condensation line, i.e. the boundary between gaseous and liquid phase of pure CO₂, and is calculated from (1). Rough limit for hydrate formation in pure water is represented by the hatched area using the data given by King [1969]. The phase boundary for hydrate formation in seawater is depicted by circles [Mitsubishi, 1990] and crosses [Ohgaki et al., 1993]. The inset shows the solubility in the pressure range 1 – 50 bar, for clarity.

0 5 10 15 20 25 0.01

0.02 0.03 0.04 0.05 0.06 0.07

~ 10 bar

~ 20 bar

~ 30 bar

~ 45.6 bar

Temperature (^oC) Cs~

calculated Stewart & Munjal (1970)

Aya et al. (1997) hydrate line

Figure 3. Solubility of CO₂ as mole fraction in seawater. The hydrate formation boundary (dashed line) is after Ohgaki et al. [1993]. Solid curves are obtained after calculating the solubility from (9) and (10) for 35 psu salinity. Open circles represent experimental data compiled by Stewart & Munjal [1970] from several different sources. Arrows show that the two outliers belong to ∼20 bar data set. Dots indicate data obtained by Aya et al. [1997] at 30 MPa in the presence of hydrate.

0 2 4 6 8 10 12 14 16 18 20

0 5 10 15 20 25 30

1 atm 10 atm

20 atm 30 atm calculated

King (1969)

Temperature (^oC)

Fugacity (atm)

Figure 4. Fugacity of CO₂ as a function of temperature and pressure. The curves are calculated using (8) and represent the fugacity of CO₂ in seawater. The circles are the values given by King [1969] for CO₂ in pure water.

2.2. Dynamics near the ocean bottom

Bottom boundary layer and turbulence

The fate of a CO2 lake disposed in the deep ocean will predominantly depend on the dynamics associated with the boundary layer near the ocean bottom and the consequent turbulence characteristics. A bottom boundary layer (BBL), a well-mixed layer and vertically uniform in potential temperature, salinity and light scattering, has been observed over large regions of the deep ocean [Weatherly and Niiler, 1974; Armi and Millard, 1976; D'Asaro, 1982]. The BBL is typically bounded by a sharp interface of strong vertical temperature, salinity and density gradients underlying a uniformly stratified region [Armi, 1977]. It was first suggested by Munk [1966], and later confirmed by Armi [1978], that mixing near the bottom may be responsible for a significant part of the vertical transport of buoyancy and tracers in the abyssal ocean. The concept has recently received a considerable amount of observational support [Thorpe et al., 1990; Lueck and Mudge, 1997; Toole et al., 1997; Ledwell et al., 2000]

and theoretical attention [Thorpe, 1987; Garrett, 1990; 1991; 2001].

Turbulence in the boundary layer is thoroughly studied [see e.g., Tennekes and Lumley, 1972; Hinze, 1975] and the oceanic BBL turbulence is widely accepted to be similar to that observed in laboratory flumes, provided that the lowermost 1-2 m part of the benthic BBL is independent of the Coriolis effects [Wimbush and Munk, 1970]. The thickness of the BBL is often taken to be about the Ekman boundary layer thickness, δE, commonly expressed by the relation

f u_*

E

= κ

δ (11)

where κ =0.41 is the von-Karman’s constant, u* is the friction velocity and f is the Coriolis parameter. For a steady turbulent boundary layer driven by the geostrophic flow, U, the friction velocity is typically estimated to be equal to U/30 [Weatherly, 1975]. Introducing U = 0.05 m s^-

1, f = 10^-4 s^-1, (11) yields δE ∼ 6.6 m, considerably less than generally observed values of BBL thickness [Armi and Millard, 1976; D'Asaro, 1982]. On dimensional grounds Weatherly and Martin [1978] proposed an expression for δE that incorporates the stratification of the fluid,

( )

[

]

E f1 N/f

A u

= +

δ (12)

where A is a constant and the buoyancy frequency in the interior, N, is given by

2 / 1

0 z

N g ⎥

⎦

⎢ ⎤

⎣

⎡

∂ ρ

∂

= ρ . (13)

They report that the BBL thickness is fairly well approximated over the range 0 ≤ N/f ≤ 200, provided that A = 1.3.

In the lower part of the BBL forces induced by shear are greater than those induced by buoyancy and Coriolis effects and the structure is similar to the ‘wall region’ of a non-rotating boundary layer where the vertical velocity distribution can be described by the logarithmic profile

0

*

z ln z ) u

z (

U = κ (14)

where z0 is a roughness parameter equal to 0.11ν/u* for smooth bottom. Roughness of the bottom can be described in terms of the roughness Reynolds number, Re* = u*d/ν, where d is the height of the roughness elements. The regime is considered to be hydraulically smooth for Re*< 5. For the present purposes, if we consider a mean flow of 0.05 m s^-1over the CO2 lake which is covered by a hydrate layer of thickness of the order of 10^-5 m (see section 4), and assuming that the height of the roughness elements is of the order of the hydrate layer

thickness, we obtain Re* ∼10^-2 using ν ∼10^-6 m² s^-1. Therefore, the flow over the CO2 lake can be considered as hydraulically smooth. Assuming the ‘law of the wall’ to be applicable –which is typically the case for smooth flow over fine sediments in lakes and oceans- the dissipation rate of turbulent kinetic energy, ε, as a function of distance from the wall, z, scales as

z u z ) U

z (

3

* 0

= κ

∂

∂ ρ

= τ

ε (15)

where the bottom shear stress, τ₀ =u_*²ρ,is used.

Vertical diffusivity

Assuming isotropy and a stationary balance between the buoyancy flux and the

dissipation of turbulent kinetic energy, ε, Osborn [1980] has proposed a scaling for the vertical diffusivity of the form

z 2

K Nε

Γ

= (16)

where Γ is the mixing efficiency and the buoyancy frequency, N, is given by (13). Γ indicates the conversion efficiency of turbulent kinetic energy into potential energy of the system. A range of estimates of Γ are found and used in the literature. Lilly et al. [1974] obtained Γ = 0.33 from atmospheric measurements. Oakey [1982] estimated Γ = 0.259 ± 0.21 using dissipation rate calculations from the high-wave number, cut-off of the temperature microstructure spectra from oceanic waters. Recently, Greenspan et al. (2001) have proposed Γ = 0.26 in a weakly stratified mixed layer. Fer et al. [2002] observed Γ = 0.15 ± 0.1 in the upper mixed layer and 0.22 ± 0.2 near the sloping sides of Lake Geneva, both comparable to the values found in the ocean interior.

Stable density stratification can strongly influence the turbulent boundary layer and for sufficiently large values of N, can suppress the turbulence completely. This has to be kept in mind, particularly with respect to strong stable density gradients induced by the dissolution of liquid CO2 in seawater, resulting in large buoyancy frequencies. Furthermore, sediment concentration, if large enough, can alter the dynamics significantly. Turbulence has to act against the gravitational forces tending to resettle the sediments which are likely to be suspended e.g. during a benthic storm.

Clauser [1956] suggested that the vertical eddy diffusion at the boundary may be estimated by

h 15u

~ 1

K_{z *} (17)

where h is the height of the bottom mixed layer, and u* is the friction velocity. Assuming that U/u* = 30, a mean velocity of U = 10 cm s^-1, and h = 50 m, (17) gives Kz ~ 10^-2 m² s^-1. This is a very rough estimate and it was pointed out by Armi [1977] that very little is known about the turbulence within the bottom boundary layer and above the Ekman layer. He suggested that roll waves, or intermittent surges, may form when the internal Froude number (ratio of the mean velocity of the current to the velocity of a long wave at the interface) of the flow is above a critical value (~1.7 for a flat bottom). Similar interfacial instabilities were observed on sloping sides of a deep lake [Fer et al., 2001; 2002]. Compared to the ocean interior, relatively large eddy diffusivities have been observed, for example Polzin et al. [1996] reported an average of Kz = 1.5x10^-2 m² s^-1 for water below 4000 m within the Romanche fracture zone and

downstream of the main sill, whereas estimated Kz in mid-depth in the vicinity was about 2x10^-

5 m² s^-1. Enhanced mixing was observed over rough topography in the abyssal ocean where measurements of tracer dispersion and turbulent energy dissipation in the Brazil Basin revealed diffusivities of 2-4 x10^-4 m² s^-1 at 500 m above abyssal hills on the flank of the Mid-Atlantic Ridge, and ∼ 10^-3 m² s^-1 nearer the bottom [Ledwell et al., 2000].

Horizontal diffusivity

Horizontal diffusivity is often written as

*

x Ahu

K = (18)

where A is a constant and h is the depth in open channel flows or the relevant vertical length scale (BBL thickness, here). Experiments in turbulent open channel flows, as well as field measurements suggest that the constant A varies greatly within the range 1 to 1000 [Fischer, 1973]. A value for Kx of O(10²) m² s^-1 can be estimated for a mean flow of 0.1 m s^-1, provided that A=600, and h = 50 m. On the Madeira Abyssal Plain, Saunders [1983] has observed

abyssal diffusivities of 250 m² s^-1 and 150 m² s^-1 for East and North directions, respectively, at a depth of 5300 m.

Elliott [1984] reported bursts of turbulent kinetic energy typically lasted for ∼50 s from turbulent current fluctuation measurements 50 cm above the seabed at a depth exceeding 4000 m on an abyssal plain. Using eddy correlation method he estimated u* = 1.2x10^-3 m s^-1 and CD = 1.9x10^-3, comparable to those used in our calculations (see section 3). D’Asaro [1982]

measured turbulence on inertial and tidal time scales which extends intermittently to the top of the boundary layer on the Hatteras Abyssal Plain.

3. The model

Temporal and spatial distribution of CO2 dissolved from a lake source in the deep ocean can be described by the unsteady, two-dimensional advection-diffusion equation which can be written in Cartesian coordinates as

( ) ( )

⎠

⎜ ⎞

⎝

⎛

∂

∂

∂ + ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

∂

∂

= ∂

∂ +∂

∂ +∂

∂

∂

z K C z x K C x z

wC x

uC t

C

z

x (19)

Here, C is the concentration of CO2, and u and w are the components of velocity in the horizontal direction, x and the vertical direction, z, respectively. Kx and Kz are the corresponding horizontal and vertical diffusivities.

We define a domain of 20 km horizontal and 200 m vertical extent located at flat ocean bottom at 3000 m depth in which there is a lake of CO2 source of 500 m length and of unit span at 1500-2000 m from the origin. The time evolution of solution to (19) is sought using finite volumes method [e.g., see Ferziger and Peric, 1999] with upwind differencing scheme for spatial and Crank-Nicolson method for temporal discretization. An equidistant grid yielding 250 m and 5 m resolution in horizontal and vertical, respectively, is employed and SIP (strongly implicit procedure) solver is used with a convergence error reduced below 10^-5. The time step is chosen to satisfy the condition of numerical stability for each of the varying cases. At the ocean floor, z = 0, a no-slip condition is imposed on the velocity, and there is no concentration flux through the ocean floor except from the source area. Experiments are terminated when the outlet boundary at x = 20 km is reached.

The hydrate film is modeled using the capillary permeation model of Mori & Mochizuki [1997]. They assume a steady state and a uniform hydrate layer with thickness, δ, in which the formation and the dissociation of the hydrate crystals are controlled by the rate of water permeation through capillaries and the rate of diffusive removal of CO2 molecules from the water side surface, respectively. Mori and Mochizuki obtained the following expression for the hydrate-film thickness, δ

( )

Samb S

2 S 2 S m

mix 2 c

C~ C~ C~

n C~ 1 nK 4

cos p

r

− +

− η

φ γ

= τ

δ (20)

where p is the porosity of the hydrate film, γ is the liquid CO2 – water interfacial tension, τ (≥1) is the tortuosity of the capillaries, φ is the water side contact angle on the capillary wall, ηmix

(~1.48x10^-3 Pa s) is the viscosity of seawater saturated with CO2, n is the hydration number and Km is the mass transfer coefficient. The solubility of CO2 in the ambient water (bulk of the water phase), C~_Samb

typically equals to zero. The thickness of the hydrate film is inversely proportional to the mass transfer coefficient, Km, which is proportional to the flow above the interface. In consequence, this suggests that δ may vary with bottom currents. Furthermore molar flux of liquid CO2 into the water,

CO2

J , can be written as

( )

⎦

⎢ ⎤

⎣

⎡

− − ρ −

= − _w

S S Samb

S mix mix m S

CO J

C~ 1

C~ C~

C~ K M

C~ 1

J ₂ 1 (21)

where ρmix and Mmix are the effective molar mass and the density of CO2 saturated seawater, respectively. The flux of water through the hydrate film, Jw can be calculated as

(

)

2 c mix mix

w C~

p 1 r M

4

J cos −

δτ ν

φ

= γ (22)

where νmix (~1.3x10^-6 m² s^-1) is the kinematic viscosity of water saturated with CO2. Mmix is defined as

S CO S

w

mix C~

M ) C~ 1 ( M

M = − + ₂ (23)

where molar mass of seawater, Mw can be estimated as 0.076 kg mol^-1 and the molar mass of CO2, M_CO2 is ∼4.4x10^-2 kg mol^-1. Mass transfer coefficient can be estimated by

Km = 0.1 u* Sc^-0.67 (24)

where Sc = ν/D is the Schmidt number (O(1000) for CO2 - seawater system).

In our calculations of the cases when a hydrate layer is present above the CO2 lake we used the parameter values of tortuosity, τ = 2, porosity, p = 10^-3,hydration number, n = 5.75, capillary radius, rc = 10^-8 m, interfacial tension, σ = 19.4x10^-3 N m^-1, and contact angle, φ = 0°.

The dependence of the hydrate film thickness on each of the parameters is depicted in Figure 5 for mean flow speeds of 0.01, 0.05 and 0.20 m s^-1. The above parameter values are retained in the calculations while the one corresponding to the relevant panel in Figure 5 is changed accordingly.

10^-7 10^-5 10^-3

1 2 3 4 5

(d) 10^-7

10^-5 10^-3

10^-1 10^-3

(c) 10^-6

10^-4 10^-2

10^-8 10^-7 10^-6

(b) 1 cm s^-1

10^-5

10^-7

0 45 90 135 180

(a) 20 cm s^-1 5 cm s^-1

Contact angle, φ ( ο)

Capillary radius, r_c (m)

Porosity, p

Tortuosity, τ

Hydrate film thickness, δ (m)

Figure 5. Dependence of the hydrate film thickness on (a) the contact angle, (b) the capillary radius, (c) the porosity, and (d) the tortuosity. The mean flow over the hydrate layer is 1 (solid line), 5 (dashed line) and 20 cm s^-1 (dotted line), respectively.

4. Results

4.1. Summary of the cases

A total number of 10 cases were simulated. They are summarized in Table 3 and the following subsections. The thickness of the BBL is set to h = 100 m for Cases 1-6, following Ohsumi et al. [1992]. This is certainly an overestimate except for the cases of benthic storm;

however, more realistic values for h are applied for the remaining cases. The vertical diffusivity above the BBL is set to 10^-5 m² s^-1, typical of the ocean interior [Toole and McDougall, 2001].

Table 3. Summary of the cases¹. Case u

(m s^-1)

Kx (m² s^-1)

Kz (m² s^-1)

hydrate density effect

1 0.05 10 0.01 absent neglected

2 0.05 10 0.01 present neglected

3 BBL² 10 BBL² absent neglected

4 BBL² 10 BBL² present neglected

5 0.20 10 0.01 absent neglected

6 0.20 10 0.01 present neglected

7 BBL³ Eq(18) BBL³ absent considered

8 BBL³ Eq(18) BBL³ present considered

9 BBL⁴ Eq(18) BBL⁴ absent considered

10 BBL⁴ Eq(18) BBL⁴ present considered

1 In all cases w = 0.

2 Boundary layer theory is used to derive the vertical distribution of velocity u and Kz for constant N. The ambient velocity above the boundary layer is chosen as 0.05 m s^-1. The mixing efficiency, Γ = 0.15.

3 Same as ² but N is calculated considering the density change effects due to dissolution of CO2.

4 Same as ³ but the ambient velocity above the boundary layer is chosen as 0.20 m s^-1. Γ = 0.25.

Case 1

The solution to (19) is sought for a constant longitudinal velocity of U = 0.05 m s^-1. It is assumed that no hydrate was formed at the interface between the liquid-CO2 lake and the seawater. The source concentration is set to 1.5 mol l^-1 which is approximately the solubility at the corresponding depth with T = 5 °C and S = 35 psu (Figure 2). Vertical and horizontal eddy diffusivities are constant and equal to 10^-2 and 10 m² s^-1, respectively. These values are

comparable to those reported in the literature (see the discussion below (17)) and agree well with (17) and (18). At each time step, the corresponding pH ( = -log[H⁺]) field is calculated for added CO2 concentration to the ambient water having T = 5 °C, S=35 psu, P = 300 bar and a constant alkalinity of 2.35 equiv. m^-3. The latter follows from assuming negligible dissolution of CaCO3 sediments. The boric acid system is included in the computations. The ambient water has a pH level of ∼7.9 at the corresponding pressure. The pH field calculation is conducted in the same fashion for all cases. The time evolution of the concentration, C, and the pH is shown in Figure 6 for 1, 2, and 3 days.

Case 2

In this case we examine the effects of the hydrate layer for the conditions of Case 1. The hydrate layer is modeled and coupled to the mean flow through (20)-(24). The values of the parameters are set to tortuosity, τ = 2, porosity, p = 10^-3,hydration number, n = 5.75, capillary radius, rc = 10^-8 m, interfacial tension, σ = 19.4x10^-3 N m^-1, and contact angle, φ = 0°. These values are retained for all the cases when the hydrate is present (i.e., Cases 2,4,6,8, and 10; see Table 3). The time evolution of the concentration, C, and the pH is shown in Figure 7 for 1, 2, and 3 days.

0 100 200

0 100 200

0 100 200

Horizontal distance (km) Height above bottom (m) ⁰

100 200

0 100 200

0 0 100 200

C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

5 10 15 20

6 7 5 4

4

5 6 7

4

5 6 7

0.1 0.05 0.2

0.05

0.05 0.1

0.1 0.2

0.2

a

b

c

d

e

f

Figure 6. Evolution of (a-c) CO2 concentration, C, distribution and the corresponding (d-f) pH

distribution for Case 1. The contours are drawn at 0.05, 0.08, 0.1, 0.15, 0.2, 0.25, 0.5, 0.75, 1, 1.25 and 1,5 mol l^-1 for C and at 1 unit intervals for pH. The CO2 lake source is located at 1.5 – 2 km from the origin.

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

Horizontal distance (km)

Height above bottom (m)

0 C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

5 10 15 20

6 7 4 5

4

5 6

7

4

5 6

7 0.1 0.05

0.05

0.05 0.1

0.1

a

b

c

d

e

f

Figure 7. Same as Figure 6, but for Case 2.

Case 3

In this case we derive a vertical distribution for the longitudinal velocity using (14) with a mean geostrophic velocity of 0.05 m s^-1 above the BBL. The thickness of the BBL is set to h = 100 m, which is about 14 times the Ekman layer thickness, δE obtained from (11) for f = O(10^-4) s^-1. A vertical profile for vertical diffusivity, Kz, is calculated using (16) with Γ = 0.15 and a constant value for N = 7x10^-4 rad s^-1, after obtaining ε profile from (15). Horizontal diffusivity is constant and the effects of density change due to the dissolution of CO2 are neglected. The time evolution of the concentration, C, and the pH is shown in Figure 8 for 1, 2, and 3 days.

Case 4

Case 4 is the same as Case 3, but the effect of the hydrate layer is included. The time evolution of the concentration, C, and the pH is shown in Figure 9 for 1, 2, and 3 days.

Case 5

In this case we examine the conditions of a typical benthic storm when no hydrate forms at the interface between the liquid CO2 and the seawater. The flow above the lake is set to 20 cm s^-1, which is a typical value for observed benthic storms (see section 5.2.). Due to large longitudinal velocity, the plume of CO2 reaches the outlet boundary after 1 d, compared to > 3 d for Cases 1-4. The time evolution of the concentration, C, and the pH is shown in Figure 10 for 12 h and 1 day.

Case 6

Case 5 is repeated with the hydrate layer present above the CO2 lake. The time evolution of the concentration, C, and the pH is shown in Figure 11 for 12 h and 1 day.

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

Horizontal distance (km)

Height above bottom (m)

0 C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

5 10 15 20

6 7 4 5

4

5 6

7

4 5 6

7 0.1 0.05

0.2

0.05

0.05 0.1

0.1 0.2

0.2

a

b

c

d

e

f

Figure 8. Same as Figure 6, but for Case 3.

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

0 100 200

Horizontal distance (km)

Height above bottom (m)

0 C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

5 10 15 20

6 7 4 5

4 5 6

7

4 5 6

7 0.05

0.1 0.05 0.05

a

b

c

d

e

f

Figure 9. Same as Figure 6, but for Case 4.

Horizontal distance (km)

Height above bottom (m)

0 5 10 15 20

7 5 6

4

5

6 7

0.05 0.05

C, 12 h

pH, 12 h C, 1 d

pH, 1 d

4

0.1

0 100 200 0 100 200 0 100 200 0 100 200

a

b

c

d

Figure 10. Evolution of (a-b) CO₂ concentration, C, distribution and the corresponding (c-d) pH distribution for Case 5. The contours are drawn at 0.05, 0.08, 0.1, 0.15, 0.2, 0.25, 0.5, 0.75, 1, 1.25 and 1,5 mol l^-1 for C and at 1 unit intervals for pH. The source of CO2 lake is located at 1.5 – 2 km from the origin.

0 100 200

0 100 200

0 100 200

0 100 200

Horizontal distance (km)

Height above bottom (m)

0 5 10 15 20

6 7 5

4

5

6 7

0.05 0.05

C, 12 h

pH, 12 h C, 1 d

pH, 1 d

4

a

b

c

d

Figure 11. Same as Figure 10 but for Case 6.

Case 7

In this case we derive a vertical distribution for the longitudinal velocity using (14) with a mean geostrophic velocity of 0.05 m s^-1 above the BBL. The thickness of the BBL is set to h = 25 m, which is about 4 times the Ekman layer thickness, δE obtained from (11) for f = O(10^-4) s^-

1. This is in agreement with oceanic observations [Armi and Millard, 1976; D'Asaro, 1982]. A vertical profile for vertical diffusivity, Kz, is calculated using (16) with Γ = 0.15, after obtaining the ε profile from (15). The density increase due to dissolution of CO2 is incorporated through a buoyancy frequency profile calculated from (13). The density of seawater with dissolved CO2, ρ2, is calculated by the relation given by Giggenbach [1990] [see also, Weiss, 1974]

) / 1000 ( C 10 x 5 . 7

C 1000

1

2 4 + ρ

= +

ρ ₋ (25)

where ρ1 is the density of seawater without dissolved CO2 [UNESCO, 1981] and C is the concentration of dissolved CO2 (here, in grams per kilogram). The strong stable stratification due to large vertical CO2 concentration gradients leads to smaller values for Kz near the bottom (see discussion of Case 8). Horizontal diffusivity is calculated using (18) with A = 10 and h = 25 m. The time evolution of the concentration, C, and the pH is shown in Figure 12 for 1, 2, and 3 days.

Case 8

Case 8 is the same as Case 7, but the effect of the hydrate layer is included. The time evolution of the concentration, C, and the pH is shown in Figure 13 for 1, 2, and 3 days.

Vertical distributions of buoyancy frequency, N, vertical diffusivity, Kz, and density, ρ, for Case 8 are shown and compared to those of Case 10 in Figure 16 both at the center of the CO2 – lake and at 10 km from the origin. Enhanced CO2 concentration close to the ocean bottom (i.e., close to the source) results in relatively high values of N, which in turn yields suppressed Kz

through (16).

Case 9

In this case we examine the conditions of a typical benthic storm when no hydrate forms at the interface between the liquid CO2 and the seawater. The velocity profile and the mixing parameters are derived as was explained for Case 7. Values used for the mean geostrophic velocity, the thickness of the BBL, and the mixing efficiency were U = 0.20 m s^-1, h = 100 m, and Γ = 0.25, respectively. The time evolution of the concentration, C, and the pH is shown in Figure 14 for 12 h and 1 day.

Case 10

Case 9 is repeated with the hydrate layer present above the CO2 lake. The time evolution of the concentration, C, and the pH is shown in Figure 15 for 12 h and 1 day.

Horizontal distance (km)

Height above bottom (m)

C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

0 5 15 10 20

0 25 50 0 25 50 0 25 50 0 25 50 0 25 50 0 25

50 a

b

c

d

e

f

5 6 7

7 7

5 6

5 6

4

4

4 0.001

0.001

0.001

Figure 12. Evolution of (a-c) CO₂ concentration, C, distribution and the corresponding (d-f) pH

distribution for Case 7. The contours are drawn at 0.001, 0.007, 0.02, 0.06, 0.2, 0.7, and 1.2 mol l^-1 (i.e.

approximately equal logarithmic space) for C and at 1 unit intervals for pH. The source of CO₂ lake is located at 1.5 – 2 km from the origin.

Horizontal distance (km)

Height above bottom (m)

C, 1 d

C, 2 d

C, 3 d

pH, 1 d

pH, 2 d

pH, 3 d

0 5 15 10 20

0 25 50 0 25 50 0 25 50 0 25 50 0 25 50 0 25

50 a

b

c

d

e

f

5 6 7

7 7

5 6

5 6

4

4

4 0.001

0.001

0.001

Figure 13. Same as Figure 12 but for Case 8

Horizontal distance (km)

Height above bottom (m)

0 5 10 15 20

C, 12 h

pH, 12 h C, 1 d

pH, 1 d

0 50 100 0 50 100 0 50 100 0 50

100 a

b

c

d

7 5

6

4

0.001

0.001 0.02

0.02 0.007

0.007

0.2 0.2

7

5

6

4

Figure 14. Evolution of (a-b) CO₂ concentration, C, distribution and the corresponding (c-d) pH

distribution for Case 9. The contours are drawn at 0.001, 0.007, 0.02, 0.06, 0.2, 0.7, and 1.2 mol l^-1 for C and at 1 unit intervals for pH. The source of CO₂ lake is located at 1.5 – 2 km from the origin.

Horizontal distance (km)

Height above bottom (m)

0 5 10 15 20

C, 12 h

pH, 12 h C, 1 d

pH, 1 d

0 50 100 0 50 100 0 50 100 0 50

100 a

b

c

d

4

0.001

4

7 6

5

7 6

5

0.001 0.007

0.007 0.02

0.02

0.2

0.2 0.06

0.06

Figure 15. Same as Figure 14 but for Case 10.

50

0.05 0

25

0.10 0.15

0

10^-1 10^-3

10^-5 Case 8 (1.75 km)

Case 8 (10.0 km) Case 10 (1.75 km) Case 10 (10.0 km)

a b

N (rad s^-1) K_z (m² s^-1)

Height above bottom (m)

1040 1045

ρ ^{(kg m-3)} c

Figure 16. Vertical distribution of (a) buoyancy frequency, N, (b) vertical diffusivity, Kz, and (c) density, ρ, for Case 8 (solid lines) and Case 10 (dotted lines). The thick lines are distributions at 1.75 km from the origin (i.e., at the center of the CO2 –lake), whereas thin lines are derived at 10 km from the origin. The profiles are obtained at 3 d for Case 8 and at 1 d for Case 10, i.e. when the CO2 plume approximately reaches the outlet boundary. The density is calculated using (25), after the CO2

distribution is obtained. Enhanced CO2 concentration close to the ocean bottom (i.e., close to the source) results in relatively high values of N, which in turn yields suppressed Kz through (16).

The total volume per unit width of water with pH < 6 is calculated and its increase with time is shown in Figure 17 for all cases. The ambient water has a pH level of ∼7.9 at the corresponding pressure.