Modeling the dissolution of settling CaCO₃ in the ocean

2.1. Model Description

[7] The water column carbonate dissolution model considers a biogenic particle comprised of either calcite or aragonite, settling through the water column. As the particle enters undersaturated waters, it experiences dissolution. A simple approach in one spatial dimension is followed. The chemistry of the bulk seawater is given as forcing data. The amount of calcite or aragonite dissolved at a certain depth is the only prognostic variable. Dissolution of biogenic carbonates in the water column depends on (1) the dissolution rate of each mineral phase, (2) the surface area in contact with water, (3) the residence time of the particles in undersaturated waters, and (4) the magnitude of undersaturation. Let C(t) (mol) be the calcite or aragonite content at time t (days). The kinetics may be formulated as [Archer, 1996]

where κ denotes the dissolution rate constant (per day) and η denotes the dissolution rate order. [CO₃²⁻] is the in situ carbonate ion concentration (μmol kg⁻¹) and [CO₃²⁻]_c is the carbonate ion concentration at saturation. An alternative rate equation for dissolution of aragonite has been proposed by Acker et al. [1987]:

where they experimentally find values of κ = 3.32 ± 1.77 d⁻¹ and η = 1.87 ± 0.15. Sensitivity studies show that application of equation (2) versus equation (1) yields almost indistinguishable results (compare Figure 3 in section 2.2). Because an analogous formulation of equation (2) for calcite does not exist, equation (1) has been used for the modeling study presented here.

[8] The carbonate ion concentration is a function of depth (i.e., pressure; the weak dependence on temperature is neglected) and is dependent on the crystal structure, namely, calcite (foraminifera and coccolithophorids) or aragonite (pteropods):

where z is depth (km), as proposed by Broecker and Takahashi [1978]. Though numerically convenient because of its simplicity, this parameterization has now been seen to be unreliable [Millero, 1995]. Thus we have fitted equations of the form of equation (3) to critical carbonate concentrations calculated by Millero [1995], yielding (see Figure 1)

[9] As a first approximation, and for use in the sensitivity experiments, bulk [CO₃²⁻] is assumed to be constant and independent of depth. This is a valid simplification, as the in situ [CO₃²⁻] is virtually constant for depths >2 km, and shallower waters are generally supersaturated with respect to calcium carbonate.

[10] Let v (m d⁻¹) be the settling velocity of the aggregate. Then, after time t the particle's depth is z(t) = vt. Thus equation (1) may be formulated as (considering calcite)

where

is the time after which undersaturation is reached and

is the dissolution rate. An analogous formulation is derived for aragonite, using the appropriate parameterization in equation (4).

[11] Standard model parameters are given in Table 1. The results of the corresponding model run are shown in Figure 2. The bulk seawater concentration of 91.7 μmol kg⁻¹ CO₃²⁻ corresponds to a calcite saturation horizon at 4 km depth and an aragonite saturation horizon at 1.65 km depth. Dissolution is defined here to be visible when 1% of the particle is dissolved. In the standard model run, dissolution of calcitic shells is visible ∼1 km below the saturation horizon. At 6 km depth, ∼68% of a settling calcite particle is dissolved within the water column. The onset of aragonite dissolution lies ∼1.3 km below the saturation horizon, as the pteropod shells sink faster than coccolithophorids. At 6 km depth, ∼57% of a settling aragonite particle is dissolved within the water column.

2.2. Sensitivity Studies

[12] This section describes results obtained with sensitivity experiments. All parameters are used with standard values except the one that is varied. Sinking rates of biogenic particles determine the time the carbonate is exposed to water column dissolution. The results of numerical experiments with settling velocities varied over 2 orders of magnitude are displayed in Figure 33. Settling velocities estimated from experimental and field work are discussed in Appendix A (Tables A1 and A2). It is noted that only for sediment depths well below the saturation horizon, significant dissolution takes place in the water column provided that the settling velocity is low. The dissolution rate is determined by the rate constant, κ, and by the order of kinetics, η, as seen in equation (1). Keir [1980] gives 3 d⁻¹ < κ < 7 d⁻¹ as a range for single coccoliths. The rate order of calcite dissolution kinetics is not well constrained, as values cited in the literature range from 1 [Hales and Emerson, 1997] to >5 [Keir, 1980] (compare Appendix B for a discussion on dissolution kinetics).

[13] The results of numerical experiments with dissolution rate constants 1 d⁻¹ ≤ κ ≤ 7 d⁻¹ are displayed in Figure 4. When compared to other parameters, varying κ has a rather minor effect on the magnitude of dissolution. Figure 5 illustrates the importance of the kinetic rate order. Varying η between 1 and 5 yields by far the largest range in model results among the sensitivity experiments. Thus better constraints on this parameter by experimental work are most desirable.

2.3. Application to Oceanic Settings

[14] Up to now, profiles with constant CO₃²⁻ bulk concentration have been used. In order to resolve the origin of the relatively shallow alkalinity maximum in the Pacific Ocean [Fiadeiro, 1980], the model is run with the in situ profile of the CO₃²⁻ concentration in this region (Figures 6 and 7). Both the Indian Ocean and the Atlantic Ocean are too supersaturated with respect to CO₃²⁻ to yield significant inorganic carbonate dissolution in the water column. The model is applied to the specific oceanic setting with standard settling velocities (10 m d⁻¹ for calcite (coccolithophorids) and 300 m d⁻¹ for aragonite (pteropods)) and with various parameterizations for the dissolution kinetics proposed in the literature (Table 2).

[15] Results for the equatorial Pacific are illustrated in Figure 8. In complete analogy the North Pacific region is considered (Figure 9). In the equatorial Pacific region, between 50 and 100% of the initial calcite is dissolved at 6 km depth, depending on the kinetic parameterization. The onset of aragonite dissolution is seen ∼0.5 km below the saturation horizon. At 6 km depth, between 55 and 96% of aragonite (pteropods) is dissolved in the water column. However, aragonite dissolution kinetics are not as thoroughly experimentally investigated as calcite dissolution kinetics, so a rate order smaller than 3 cannot be excluded, which would lead to increased amounts of aragonite dissolution in the water column.

[16] The North Pacific model run shows increased dissolution mainly for calcite: Dissolution starts close to the saturation horizon at ∼1.2 km water depth. Between 2 and 3 km depth, dissolution slows down due to the near-equilibrium situation (Figure 7). With linear kinetics, nearly all calcite is dissolved at 4 km depth, which is at the onset of calcite dissolution in the equatorial Pacific. Aragonite dissolution reaches 100%. The parameterizations of Morse et al. [1979] and of Byrne et al. [1984] yield almost indistinguishable results.

[17] The drastic difference in calcite dissolution between the equatorial and North Pacific is due to the location of the calcite saturation horizon. As can be seen in Figures 6 and 7, the saturation horizons of calcite and aragonite are separated by ∼3 km in the equatorial Pacific, while they are both shallower than 1 km in the North Pacific.

3.1. Model Description

[19] The aggregate is modeled as an idealized porous sphere consisting of water, organic carbon, and calcite. Initial data are the C_org content (μg C), the calcite content (μg CaCO₃), and the porosity. The volume of the aggregate is then calculated as

where ϕ is the porosity and V_c and V_cal (m³) are the volumes of the organic carbon and calcite content, respectively, which are derived from mass via the specific densities. The standard value for porosity is ϕ = 0.9, following in situ observations by Alldredge and Gotschalk [1988]. The settling velocity v (m d⁻¹) is calculated via the parameterization by Alldredge and Gotschalk [1988], who observed settling velocities of marine snow aggregates to be less than those of an equivalent spherical particle settling according to Stokes' law:

where d_a = 2r_a (mm) is the diameter of the aggregate.

[20] The CO₂ produced by organic matter respiration is assumed to be released at the spherical surface of the aggregate. Dissolution via chemical alteration of the microenvironment requires the establishment of a stable diffusive boundary layer around the sinking particle, for which the aggregate radius and the sinking velocity are crucial. The concentration C_a of a dissolved species (i.e., CO₂, CO₃²⁻, NO₃⁻, or O₂) at the aggregate surface may be calculated using the relation [Ploug et al., 1999b]

where Q (mol s⁻¹) is the area-integrated flux of the dissolved species in question out of the aggregate and D (m² s⁻¹) is the appropriate diffusion coefficient. C_∞ and C_a (mol m⁻³) are the concentration of the species in the bulk medium and at the aggregate surface, respectively, and δ_eff (m) is the effective thickness of the diffusive boundary layer. Because of shear, the latter decreases with increasing settling velocity. The relative increase of advection compared to diffusivity is described by the Sherwood number Sh [Sherwood et al., 1975]:

In a stagnant fluid, Sh = 1, and thus δ_eff = r_a. The Sherwood number varies with both the Reynolds number Re and the Péclet number Pe:

where ν is the kinematic viscosity of sea water (∼0.7 × 10⁻⁶ m² s⁻¹).

[21] Parameterizations of the Sherwood number are given by Karp-Boss et al. [1996]:

[22] Solving equation (10) for C_a and using equation (11) yields

Applying equation (14), we calculate ΣCO₂ and alkalinity concentrations at the aggregate surface, considering the fluxes deriving from remineralization of organic matter and dissolution of carbonate. From ΣCO₂ and alkalinity, [CO₃²⁻] is derived via the carbonate dissociation constants from Department of Energy [1994] and using pressure corrections by Culberson and Pytkowicz [1968].

[23] In the case of undersaturation, calcite is dissolved as described in section 2.1 (equation (1) and Table 1). Organic carbon is remineralized according to first-order kinetics,

where λ (d⁻¹) is the remineralization rate constant; that is, 1/λ gives the e-folding time. Respiration of 1 mol organic carbon consumes moles oxygen, where is the Redfield ratio of oxygen to carbon (compare Table 3). Using equation (14) and bulk oxygen concentrations from GEOSECS data [Takahashi et al., 1981], the maximal amount of C_org is calculated that may be oxidized in one time step, avoiding oxygen depletion.

Table 3. Parameter Values

Symbol	Parameter	Value	Reference
ϕ	porosity	0.9	Alldredge and Gotschalk [1988]
	Corg content	0–15 μg C	compare Table 4
	calcite content	0–3 μg	Beers et al. [1986]; Tyrell and Taylor [1996]
	density calcite	2.7 g cm−3	Young [1994]
	density Corg	1.07 g cm−3	Young [1994]
κ	calcite dissolution rate constant	5 d−1	Keir [1980]
η	calcite dissolution rate order	4.5	Keir [1980]
λ	Corg decomposition time constant	1/15 d−1	Westrich and Berner [1984]
	Redfield ratio O2 to C	170/117	Anderson and Sarmiento [1994]

[24] The numerical implementation of the model is as follows. Initial conditions are the content of calcite and organic carbon. The radius of the aggregate and initial settling velocity are computed. The time step is set to 6 hours. At each time step, for each dissolved species, the Sherwood number is computed from the radius and settling velocity. Then the fluxes of CO₂ and NO₃⁻ are derived from the remineralization of organic matter. In the case of undersaturation the flux of CO₃²⁻ is calculated. From these fluxes and the bulk chemistry, ΣCO₂ and TA at the aggregate surface are derived according to equation (14), setting the new value for calcite saturation. The radius of the aggregate is decreased according to the amounts of organic matter respired and calcite dissolved, while keeping porosity constant. Then the new settling velocity and depth are computed. This scheme is iterated until the aggregate has reached the seafloor at a prescribed depth.

3.2. Parameter Values

[25] The model parameter values are given in Table 3. Experimentally derived parameterizations of the calcite dissolution kinetic are discussed in Appendix B. Westrich and Berner [1984] give a compilation of experimentally determined first-order decay constants for oxic organic matter decomposition. The e-folding time varies between 11 and 122 days, with a mean of 32 days. Walsh et al. [1988] find an average e-folding time of 38 days in sediment traps deployed in three North Pacific stations. Experimentally, Westrich and Berner [1984] derive e-folding times of 15 days for oxic decomposition; this value has been used for the model runs.

[26] The aggregate's calcite content is considered to be Emiliania huxleyi cells. In a field study of two E. huxleyi blooms in the Gulf of Maine, Balch et al. [1991] find 25–100 and 100–400 coccoliths per cell. Fernández et al. [1993] point out that numbers are much lower in open ocean blooms, as they find 20–40 coccoliths per cell in an E. huxleyi bloom in the northeast Atlantic Ocean. Tyrrell and Taylor [1996] state that 10–50 coccoliths per cell is typical, with one coccolith containing 1.7–8.3 pg calcite. Trent [1985] and Beers et al. [1986] observed coccolithophorid cell numbers per aggregate ranging between 180 and 870. Assuming 30 coccoliths per cell, 5 pg calcite per coccolith, and 1000 cells per aggregate yields 0.15 μg calcite (≅1.8 × 10⁻² μg C). There are hints that at the end of the bloom, when marine snow aggregates are formed, the coccolithophorids are more heavily plated [Balch et al., 1991]. Adopting 100 coccoliths per cell and 5000 cells per aggregate gives 2.5 μg calcite content at the upper end of the range.

[27] The organic carbon content of marine snow aggregates generally ranges from 0.1 to 12 μg C, but may be as high as 220 μg C (Table 4). The calculations presented here adopt values of 1–15 μg C per aggregate. It has to be stressed that there is no functional relationship between organic carbon content and the number of cells of calcareous organisms contained in a marine snow particle.

3.3. Model Runs

[28] In order to investigate the capacity of respiration-driven dissolution, a number of model runs are performed (Table 5). Figures 10 and 11 show results of a model run with GEOSECS data of the North Pacific [Takahashi et al., 1981] determining the bulk carbonate system. A marine snow aggregate is considered, consisting of 15 μg C_org and 1000 E. huxleyi cells, with 0.15 μg calcite.

[29] The initial diameter (Figure 11b) is 1.75 mm, which lies in the range of marine snow aggregates observed by Alldredge and Gotschalk [1988]. The settling velocity is initially 58 m d⁻¹ and decreases down to 42 m d⁻¹ (Figure 11a). The Sherwood numbers decrease from 4–5 down to 2–3. This results in the diffusive boundary layer thickness ranging between 25 and 200 μm, which conforms to work by Alldredge and Cohen [1987], who measured chemical boundary layers hundreds of micrometers thick in marine snow aggregates in a flow tank, simulating advective forces as they appear when the aggregate is settling through the water column. Respiration of the organic matter increases [CO₂] at the aggregate surface over the bulk concentration by a maximum of 7 μm kg⁻¹ (Figure 10a), which lowers the calcite saturation, as is visible in Figure 10c. However, because of rather fast settling velocities and slow dissolution kinetics, this does not have a significant effect on calcite dissolution. At a depth of 3 km, nearly all organic matter has been respired, so the amount of calcite dissolution at this depth is a good measure for dissolution attributable to respiration effects. Here, a mere 0.1% of the calcite content is dissolved, in contrast to no dissolution in the absence of organic matter remineralization.

[30] The model is then tested for different dissolution kinetics proposed in the literature (Table 2). As Keir's [1980] is the only work specifically dealing with coccoliths, it is assumed to be the most adequate kinetic for the present model. However, other results should not be excluded; for example, the work of Walter and Morse [1985] deals with low-Mg calcite, which is attributable to coccoliths. At 3 km depth, no calcite dissolution is evident in the equatorial Pacific, no matter what kinetic parameters are chosen. Figure 12 shows the same model run as Figure 10, but with dissolution kinetics proposed by Morse [1978] (κ = 5 d⁻¹, η = 3). In this case, 4% of the initial calcite content is dissolved at 3 km depth, compared to 0.04% without organic matter remineralization.

[31] Table 6 displays results varying the amount of organic matter in the aggregate. Inclusion of organic carbon always enhances the amount of calcite dissolved. In case C, where κ = 5 d⁻¹ and η = 3, the influence of respiration is visible. In the case of linear kinetics (case D) only, where κ = 0.38, d⁻¹, and η = 1, calcite dissolution due to respiration processes is significant.

[32] Table 7 shows model results for the same model run as in Table 6, but with 2.5 μg calcite content. Note that with a given dissolution rate a certain fraction rather than a certain amount of the calcite present is dissolved over some time interval (compare equation (1)). In principal, when doubling the amount of calcite present, mass dissolution is expected to double, too. Only when the release of carbonate ions from dissolution significantly perturbs the carbonate system toward a more saturated state is the amount of calcite present relevant for the fraction dissolving per time interval. Therefore incorporating more calcite into the aggregate should not lead to less dissolution fraction wise simply because there is more calcite present, but because more calcite should make the aggregate larger and should make it settle faster.

[33] However, the fraction dissolved at a depth of 3 km does not differ significantly from the model results displayed in Table 6. This is due to the fact that calcite, with its large density (2.7 g cm⁻³), does not influence the radius and hence does not influence the settling velocity of the aggregate significantly. This might be a drawback of the model formulation. Nevertheless, other estimates on marine snow settling velocities (see Table A2) conform to the settling velocities calculated by the present model, though they are a bit higher, as most aggregates investigated are between 1 and 10 mm in diameter, and the aggregate sizes modeled here are on the lower end of this range. The marine snow particles investigated by Alldredge and Gotschalk [1988] to derive the formulation (equation (9)) partly contained diatoms or fecal pellets, which are comparable in density to calcite.

[34] We now investigate what parameter variations would yield an aggregate that exhibits significant calcite dissolution. The main reason for the low effectiveness of respiration-driven dissolution is the porosity of the aggregate. High porosities lead to large diameters. Given a constant flux of, for example, CO₂, the concentration gradient between the aggregate surface and the bulk medium decreases with increasing radius. Thus larger diameters at constant source strength Q are unfavorable for creating locally undersaturated chemical environments. To overcome these disadvantages with respect to calcite dissolution, the aggregate has to be loaded with large amounts of organic matter, which in turn increases the radius drastically due to the low density of C_org. This makes the aggregate settle faster, according to equation (9). The effect of a high settling velocity is twofold. It reduces the time the aggregate spends in undersaturation, and it decreases the thickness of the boundary layer, which in turn decreases the concentration gradient, as can be seen from equations (11)–(14).

[35] With small radius and resulting low settling velocity, the thickness of the diffusive boundary layer is approximately r_a, while it is < r_a at large radius and high settling velocities due to the increase of the Sherwood number with increasing velocities. This promotes loss of respiratory CO₂ and thus less dissolution. Because of the high porosity of the aggregate and its resulting fractal geometry, the mass, i.e., organic carbon content, does not scale as r_a³, but in fact only as ∼r^1.4, corresponding to the fractal dimension [Logan and Wilkinson, 1990; Alldredge, 1998]. Thus the CO₂ source cannot keep up with the increasing surface (scaling as r_a²) through which it diffuses when the aggregate gets larger. Thus, given the high porosity of the aggregate, an increase in size is unfavorable for respiration-driven dissolution because of increased settling velocity (i.e., less time spent in the water column), diminished thickness of the boundary layer (resulting in decreased concentration gradients), and a negative scaling ratio of respiratory CO₂ source to diffusion area.

[36] Because of the idealized geometry of the aggregate in the model, the organic carbon content in the model actually scales as r_a³. Nevertheless, respiration-driven dissolution correlates negatively with size, which means that the other two mechanisms outweigh the (in the model) positive scaling ratio of respiratory CO₂ source to diffusion area. Thus the model somewhat overestimates the potential of respiration-driven dissolution of CaCO₃ in marine snow aggregates.

[37] In order to quantify these relationships, the model is run with variable porosity. Results are illustrated in Figure 13. Respiration-driven dissolution is significant in the North Pacific provided ϕ ≤ 0.4 and is significant in the equatorial Pacific when ϕ ≤ 0.2. However, these results may be over estimations, since the settling velocity parameterization (equation (9)) has been derived with aggregate shaving porosity ϕ ≥ 0.9. Aggregates with smaller porosity are expected to settle faster, because the density increases with decreasing porosity. Stokes' law is not valid for the model's parameter range in Figure 13 since its application leads to Reynolds numbers >1.

[38] A specific model run with ϕ = 0.3 and 10 μg C organic carbon is shown in Figure 14. Between 1 and 2 km depth, 14% of the calcite content is dissolved due to organic matter respiration. Then the organic matter has vanished, and calcite dissolution ceases until, at ∼4.5 km depth, bulk undersaturation is strong enough to trigger some more dissolution. The gradients in [CO₂] and [CO₃²⁻] (Figures 14a and 14b) are much larger than with porosity 0.9 (compare Figure 10). The kink at 1 km is due to reduced organic matter remineralization as a consequence of oxygen deficiency, which is illustrated in Figure 15. Decreasing the porosity reduces the diffusive fluxes and the aggregate diameter, which enhances the oxygen concentration gradient. Oxygen deficiency appears with 5, 10, and 15 μg C_org when ϕ ≤ 0.3, 0.5, and 0.625, respectively.