Supporting Information for
“Catalysis and Chemical Mechanisms of Calcite Dissolution in Seawater”
Adam V. Subhas, Jess F. Adkins, Nick E. Rollins, John Naviaux, Jonathan Erez, and
William M. Berelson
Calculating a Reactive Layer Thickness using Secondary Ion Mass
Spectrometry (SIMS)
Data Collection
6-8 large (300x300
μ
m)
13
C-labeled calcite grains were embedded in 1-inch epoxy discs. Discs
were then polished to 0.25
μ
m and reacted for 48 hours in an experimental setup similar to
the dissolution experiments described in the Methods section and in Subhas et al. (2015). The
alkalinity contribution of the epoxy resin to alkalinity over the course of the 48-hour experiment
was minimal (at most a 4
μ
eq/kg increase). Omegas for the supersaturated and undersaturated
experiments were 1.32
±
0.03 and 0.96
±
0.02. After reaction, discs were rinsed thoroughly with
methanol and air-dried for 30 minutes before covering and storing in the dark. Methanol was
used to avoid further reaction with the surface of the calcites, and for its low vapor pressure
(thus low potential for carbon contamination during SIMS analyses).
Secondary Ion Mass Spectrometer profiles were collected on a Cameca 7f SIMS instrument.
Carbon isotope ratios were collected in positive ion mode using a Cesium beam. Mg/Ca ratios
were collected as
24
Mg
/
40
Ca in negative ion mode using an oxygen beam. A correction factor
of 1.23 was applied to collected Mg/Ca ratios to account for the relative natural abundances
of these isotopes. Isotope fractionation factors for both carbon and Mg/Ca incorporation into
calcite were ignored as they are on the order of a permil relative and thus insignificant at our level
of precision. Pertinent tuning parameters can be found in Table S1. Before measuring isotopic
profiles, the epoxy discs were sputter-coated with 20 nm gold to ensure surface conductivity.
Coated discs were then loaded into the sample evacuation chamber for
∼
30 hours to degas the
epoxy. Dynamic transfer optics setting (DTOS) was turned off to minimize edge effects on the
profile ion beam.
Profiles were collected on at least 3 different grains, with at least 2 different spots on each
grain for a total of at least 6 total analyses per experiment. Profiles on unreacted grains were
collected in the same way without reaction in seawater. SIMS profiles in the main text show
carbon isotope profiles for each epoxy disc. Solid lines are the mean of all profiles, and shaded
regions depict the standard deviation of all profiles collected under each experimental condition.
Proflies for Mg/Ca are shown in Figure S1; profiles of carbon isotope ratios are shown in Fig.1
in the main text. Carbon isotope profiles are generally cleaner and show greater separation
between supersaturated and undersaturated experiments. Mg/Ca profiles are noisier, and there
is less distinction between the supersaturated and undersaturated experiments. This difference
can be explained by two mechanisms. First, fewer counts of Mg than Carbon-13 led to a greater
relative error in each cycle, calculated using Poisson statistics. Especially deeper in the profile,
this counting error manifests as noisy profiles. Second, cesium and oxygen ion beams have
different properties and thus mixing dynamics during the sputtering process. The oxygen beam
1
was more energetic and produced rougher pits compared to the cesium beam (Figure S2).
Plots of SIMS profiles were corrected for the thickness of gold, which was burned through in
about 12 cycles for C isotopes and 7 cycles for Mg/Ca. The appearance of CaCO
3
underlying
the gold coating was diagnosed by an increase in the raw counts of C and Ca ions. In order to
convert mass spectrometer cycles into vertical depth, the pit depths of all profiles were measured
using a profilometer, and applied as a mean value to all profiles. This mean profile depth was
used to convert the number of SIMS analysis cycles to depth, corrected for the thickness of
the gold coating. Examples of pit geometries are shown in Figure S2. Pits were in general
shallower and rougher in topography for the oxygen beam (Mg/Ca) than the cesium beam (C
isotope) profiles. This difference in profile shape helps to explain the different profile shapes
seen between Figure 1 in the main text and Figure S1. Poorer resolution in the depth direction,
and a more energetic oxygen beam, lead to poorer distinction between magnesium signatures in
the supersaturated and undersaturated experiments in Figure S1.
Analysis and Data Reduction
Because SIMS sputtering is an energetic process, profiles shown in Figure 1 and S1 represent
mixing curves, where excavated calcite is continually mixed down into the solid during the
sputtering process. The number of moles in the reactive calcite layer was estimated using both
cation and anion mass balances by integrating the amount of
12
C and Mg throughout the profile.
Both carbon isotope and Mg/Ca intensity ratios were converted to intensity fractions for each
ratio cycle, and multiplied by the number of moles in each cycle to give a number of either Mg
or
12
C moles in each cycle. These values were then summed to give a total number of Mg
2+
and
12
C moles added in the solid:
12
C
tot
=
i
∑
0
z
i
·
(
SA
raster
·
ρ
calcite
mm
calcite
)
·
f
12
i
;
Mg
2+
tot
=
i
∑
0
z
i
·
(
SA
raster
·
ρ
calcite
mm
calcite
)
·
f
Mg
2+
i
;
(1)
where the sum was taken over all cycles
i
in the profiles.
z
i
is the mean thickness of calcite
sputtered in each cycle, estimated by dividing the measured pit depth by the total number of
cycles.
f
12
i
is the intensity fraction of
12
C
/
(
13
C +
12
C) measured in SIMS cycle
i
;
f
Mg
2+
i
is
the intensity fraction of Mg
2+
measured in each SIMS cycle.
SA
raster
is the SIMS raster area
(50x50
μ
m
2
);
ρ
calcite
= 2
.
6
·
10
6
g m
−
3
is the density of calcite and
mm
calcite
= 100 g mol
−
1
is
the molar mass of calcite. For carbon isotope balance, any
12
C enrichment above the control run
was assumed to represent new calcite added. Thus, the mean integrated number of moles from
the control experiment was subtracted from the experimental profiles. For Mg/Ca balance, an
Mg/Ca ratio of 0.08 was assumed for the newly precipitated solid (Mucci and Morse, 1984; Sun
et al., 2015). Thus to calculate the number of moles of new calcite added we divided by this
ratio:
2
CaCO
3
(
12
C
,new
)
=
12
C
tot
;
CaCO
3
(
Mg,new
)
= Mg
2+
tot
·
Ca
Mg
∣
∣
∣
∣
solid,eq
.
(2)
Here, Ca
/
Mg
solid,eq
= 1
/
0
.
08 is taken from measurements of about 8 mole % magnesium
incorporated into calcite grown in equilibrium with modern seawater (Mucci and Morse, 1984).
Estimates for the total moles of new/exchanged calcite, as calculated by carbon and magnesium
mass balances, as well as effective reactive layer thicknesses, are presented in LayerThickness.
Effective layer thicknesses (
z
ef f
) were estimated by converting the total moles of new calcite to
a volume of calcite (
ρ
calcite
= 2
.
63
gcm
−
3
) and dividing by the sims raster surface area (50x50
μ
m
2
):
z
ef f,
12
C
= CaCO
3
(
12
C
,
new)
·
mm
calcite
ρ
calcite
·
SA
raster
;
z
ef f,Mg
= CaCO
3
(
Mg,new
)
·
mm
calcite
ρ
calcite
·
SA
raster
.
(3)
The values for
z
ef f,Mg
are lower than
z
ef f,
12
C
. We do not know the precipitated phase in our
system, and thus there is significant uncertainty in the newly precipitated solid Mg/Ca ratio.
For example, there is experimental evidence that in seawater of Mg/Ca = 5, aragonite is favored
to precipitate even on calcite seeds (Sun et al., 2015; Morse et al., 1997). The distribution
coefficient of Mg
2+
into aragonite is significantly lower than that in calcite by at least one order
of magnitude (Zhong and Mucci, 1989). An Mg/Ca of 0.08 is thus an upper bound on the
amount of Mg incorporated into our calcite. Using it in Eq.(2) represents a lower limit on the
total carbonate precipitation in these experiments, and helps explain why layer thicknesses using
Mg/Ca mass balance are thinner than those from carbon isotope mass balance.
Finally, we ground-truthed the SIMS data using the precipitation (Ω = 1
.
3) experiment
and compared to literature values of precipitation rate. At Ω = 1
.
3, Zhong and Mucci (1989)
measured a calcite precipitation rate of 1
μ
mol m
−
2
hr
−
1
. Mulitplying this rate by the molar
mass and dividing by the density of calcite gives a rate of 0.037 nanometers of calcite precipitated
per hour. Over a 48 hour period, we should thus expect about 1.8 nm of calcite precipitated
onto our SIMS disks. However, we measured that about 3-9 nm of
12
C-calcite was added in our
Ω = 1
.
3 experiments (Table S2). The discrepancy between the net amount of calcite added to
the surface and the amount of
12
C measured on the SIMS can be attributed to calcite exchange
due to dissolution-precipitation reactions at the calcite surface.
Assuming a linear relationship between precipitation rate and Ω, similar to Eq. 1 in the main
text, there should be a 30% imbalance in gross precipitation and dissolution rates at Ω = 1
.
3.
Dissolution and reprecipitation will exchange solid
13
C with seawater
12
C, above and beyond
the net accumulation of calcite. To estimate the amount of
12
C added in such dissolution-
precipitation reactions, we used the box model described below, but modified it to add calcite
to the reactive surface at a rate of 1
μ
mol m
−
2
hr
−
1
. We also set the balance of dissolution to
precipitation of
R
precip
/R
diss
= 1
.
3. The results of this model, and a comparison to the model
run in its “dissolution” configuration, are shown in Figure S3. As seen in Figure S3b, about 2.1
nanometers of calcite was added under precipitation conditions. Figure S3a further shows that
about 4.4 nanometers of
12
C-calcite was added to the solid through a combination of dissolution
3
and precipitation reactions. 4.4 nanometers is within our integrated
12
C measurement of 3.3-9.1
nanometers of
12
C-calcite added in SIMS experiments at Ω = 1
.
3. The model thus confirms our
SIMS measurements are accurately measuring the amount of precipitation in the solid, and gives
us confidence in our measurement of the
12
C enrichment in our undersaturated experiments.
Dissolution-Precipitation Box Model and Data Fitting
Box Model
We developed a box model of dissolution and precipitation at the mineral surface to model our
raw data of
13
C versus time. A schematic of this model is shown in Figure S4. The model has
three main reservoirs: The bulk solution, a boundary layer, and a reactive calcite layer. Carbon
from the bulk solution diffuses into and out of the boundary layer. Dissolution and precipitation
reactions occur between the boundary layer and the reactive calcite layer. Because so little
calcite is dissolved in our experiments, we keep the total surface area reacting with seawater
fixed over time. We model this reactive calcite layer as a fixed volume “reaction front” that is
constantly in contact with seawater, calculated as the fixed total surface area multiplied by a
reactive layer thickness. In order to keep this volume fixed through time in the model, material
removed by dissolution must be replaced by an equal flux of material from the calcite interior
into the reaction front:
d
[
13
C]
bulk
dt
=
D
SA
l
(
[
13
C]
bl
−
[
13
C]
bulk
)
1
V
bulk
;
(4a)
d
[
12
C]
bulk
dt
=
D
SA
l
(
[
12
C]
bl
−
[
12
C]
bulk
)
1
V
bulk
;
(4b)
d
[
13
C]
bl
dt
=
(
D
SA
l
(
[
13
C]
bulk
−
[
13
C]
bl
)
+
R
diss
·
f
13
solid
−
R
precip
·
f
13
bl
)
1
V
bl
;
(4c)
d
[
12
C]
bl
dt
=
(
D
SA
l
(
[
12
C]
bulk
−
[
12
C]
bl
)
+
R
diss
·
f
12
solid
−
R
precip
·
f
12
bl
)
1
V
bl
;
(4d)
d
[
13
C]
solid
dt
=
R
precip
·
f
13
bl
−
R
diss
·
f
13
solid
+ (
R
diss
−
R
precip
)
·
f
13
interior
;
(4e)
d
[
12
C]
solid
dt
=
R
precip
·
f
12
bl
−
R
diss
·
f
12
solid
+ (
R
diss
−
R
precip
)
·
f
12
interior
.
(4f)
Fluxes in solution (Eqs. 4a-4d) were calculated in units of concentration per time (moles
m
−
3
sec
−
1
). The subscripts
bulk
and
bl
correspond to the bulk solution and diffusive boundary
layer reservoirs, respectively. Fluxes into and out of the solid (Eqs. 4e-4f) were calculated in
units of moles per time.
f
i
represent the isotopic mole fractions of stable carbon isotope
i
. Rates
of precipitation (
R
precip
) and dissolution (
R
diss
) are in units of moles per time, and represent the
total amount of precipitation or dissolution. The third terms in Eqs. 4e and 4f model the supply
of new calcite from the interior, at a rate equal to the net dissolution rate (
R
diss
−
R
precip
). This
material carries an isotopic composition of the calcite interior. In this case our calcite is pure
13
C, so
f
13
interior
= 1 and
f
12
interior
= 0.
4
The following mass balance constraints were applied to the above differential equations:
N
solid,total
=
m
calcite
·
SA
calcite
·
z
mono
·
n
mono
·
ρ
calcite
mm
calcite
;
(5a)
f
13
solid
=
N
13
solid
N
solid,total
;
(5b)
f
12
solid
= 1
−
f
13
solid
;
(5c)
V
bl
=
SA
calcite
·
z
bl
;
(5d)
f
13
bl
=
[
13
C]
bl
[
13
C]
bl
+ [
12
C]
bl
;
(5e)
f
12
bl
= 1
−
f
13
bl
,
(5f)
where
N
solid,total
is the total number of moles in the reactive calcite layer;
V
bl
is the volume in
cubic meters of the boundary layer and is equal to the calcite surface area
SA
calcite
multiplied
by the boundary layer thickness
l
, assumed to be 10
μ
m.
SA
calcite
= 0
.
09 m
2
g
−
1
was measured
using Kr-BET (Subhas et al., 2015). The thickness of a monolayer of calcite was assumed to be
z
mono
=0.5 nm, and the number of monolayers
n
mono
was varied as discussed below to change
the overall volume of calcite reacting with solution.
mm
calcite
= 100 g mol
−
1
is the molar mass
of calcite and
ρ
calcite
= 2
.
63 g cm
−
3
is its density.
The bulk volume used is identical to our experimental conditions (300 grams of seawater at
density 1025 kg m
−
3
)
.
R
precip
=
R
diss
/r
f b
, the ratio of dissolution to precipitation, which varied
between 1.001 to 10 or more (see below for model-data fitting). The model was developed and
run in MATLAB r2015a. If the boundary layer
l
is very thick, diffusion out of the boundary
layer restricts the expression of the curvature in the bulk data (not shown), because the initial
burst of
13
C-labeled DIC is slowly released and mixed out of the boundary layer, rather than
being delivered to the bulk solution immediately. However, we expect our boundary layer to be
relatively thin and do not expect diffusion to be a major component of our observed rates. This
is because mixing rate does not affect our measured dissolution rates between 60-90 rpm. A
10
μ
m boundary layer and the carbonate ion diffusion coefficient in seawater (9.55
·
10
−
10
m
2
s
−
1
,
Li and Gregory (1974)) expresses the curvature in
δ
13
C versus time in both the boundary layer
and the bulk solution. Example model output data, taken from the same model run as Figure 2a
in the main text, is shown in Figure S5 for the bulk solution, boundary layer, and solid reactive
calcite reservoirs. The curvature in all three plots corresponds to the solid calcite layer coming
into steady state with respect to the dissolution and precipitation fluxes.
A mass balance can also be constructed for the Mg/Ca of calcite in solution using similar
5
equations to Eq.(4):
d
[Mg
2+
]
bulk
dt
=
D
SA
l
(
[Mg
2+
]
bl
−
[Mg
2+
]
bulk
)
1
V
bulk
;
(6a)
d
[Ca
2+
]
bulk
dt
=
D
SA
l
(
[Ca
2+
]
bl
−
[Ca
2+
]
bulk
)
1
V
bulk
;
(6b)
d
[Mg
2+
]
bl
dt
=
(
D
SA
l
(
[Mg
2+
]
bulk
−
[Mg
2+
]
bl
)
+
R
diss
·
f
Mg
2+
solid
−
R
precip
·
f
Mg
2+
bl
)
1
V
bl
;
(6c)
d
[Ca
2+
]
bl
dt
=
(
D
SA
l
(
[Ca
2+
]
bulk
−
[Ca
2+
]
bl
)
+
R
diss
·
f
Ca
2+
solid
−
R
precip
·
f
Ca
2+
bl
)
1
V
bl
;
(6d)
d
[Mg
2+
]
solid
dt
=
R
precip
·
f
Mg
2+
bl
−
R
diss
·
f
Mg
2+
solid
+ (
R
diss
−
R
precip
)
·
f
Mg
2+
interior
,
(6e)
d
[Ca
2+
]
solid
dt
=
R
precip
·
f
Ca
2+
bl
−
R
diss
·
f
Ca
2+
solid
+ (
R
diss
−
R
precip
)
·
f
Ca
2+
interior
,
(6f)
Where
f
Me
2+
i
is the mole fraction of either Mg
2+
or Ca
2+
in the reservoir
i
. As in Eq.(4),
solution fluxes are in terms of concentration (moles per volume), and solid fluxes are in terms of
total moles.
R
diss
and
R
precip
are the total rates of dissolution and precipitation, respectively,
in units of moles per time. Incorporation of Mg
2+
into the solid is sensitive to the ratio of
magnesium to calcium in solution, and the distribution coefficient of magnesium into calcite.
We modified
f
Mg
2+
bl
such that it represents the mole fraction of magnesium precipitated from
the solution:
f
Mg
2+
bl
=
D
Mg
·
(Mg
/
Ca)
bl
1 +
D
Mg
·
(Mg
/
Ca)
bl
.
(7)
Here, D
Mg
= (Mg
/
Ca)
solid
/
(Mg
/
Ca)
bulksolution
= 0
.
019 (Oomori et al., 1987) is the distri-
bution coefficient of Mg
2+
into calcite at room temperature, and (Mg
/
Ca)
bl
is the Mg/Ca ratio
in the model boundary layer. The quantity
D
Mg
·
(Mg
/
Ca)
bl
is thus the Mg/Ca ratio of the
precipitated solid in equilibrium with calcite. The right hand side of this equation converts
mole ratio to mole fraction, analogous to converting carbon isotope ratios into carbon isotope
mole fractions. Implicit in this equation is the assumption that at every model time step, cal-
cite is precipitating in equilibrium with the solution. Implementation of cation mass balance
in the model is also shown in Figure S5. As expected, the Mg/Ca ratio of the solid decreases
as dissolution overtakes precipitation. At equilibrium, the solid Mg/Ca =
D
Mg
·
(Mg
/
Ca)
bulk
∼
0
.
09
.
Data Fitting Using the Box Model Output
Data fitting was performed using the grid search method. A large model parameter space was
generated, and data were compared to all model runs over the entire parameter range. The
model was initialized and run for 400 values of
R
diss
and 400 values of
r
f b
, spanning multiple
orders of magnitude of dissolution rate (10
−
13
to 10
−
8
moles/s), and varying
r
f b
from 1.001-10 (in
the framework above, 1-Ω from
∼
0.009 to 0.9). The time-evolution of each of these model runs
(i.e. moles dissolved versus time) was then stored. Experimental data were screened to make
sure that there was sufficient data density to provide good constraints on dataset curvature. For
6
each experiment, the initial data point was set to zero moles dissolved, to eliminate excess
δ
13
C
increase that may have occured due to fine particle dissolution. Model data were interpolated
to the time values of the experimental data points, and the normalized fitting parameter F was
calculated as:
nF
=
m
net
−
m
m,net
m
net
+
1
i
i
∑
t
=0
|
n
t
−
n
m,t
|
n
t
;
(8)
where in the first term,
m
net
is the experimentally determined net rate, and
m
m,net
is the
modeled net rate, calculated as:
m
m,net
=
R
diss
−
R
precip
=
R
diss
−
R
diss
r
f b
=
R
diss
(1
−
1
r
f b
)
,
(9)
where
R
diss
and
r
f b
are defined as above. In the second term,
n
t
are the measured number of
moles dissolved at time points 0 to
i
, and
n
m,t
are the model-calculated moles dissolved at the
interpolated time points 0 to
i
. A normalized fitting parameter was used to avoid bias in the
fitting routine arising from differences in the absolute magnitude of dissolution. The implicit
assumption of these two terms is that they are both weighted equally in calculating the fitting
parameter nF. The fitting parameter was calculated for each dataset at all 400x400 values of
R
diss
and
r
f b
. The best-fitting values of
R
diss
and
r
f b
to each experiment dataset were then
found by identifying the
R
diss
and
r
f b
values which gave the lowest fitting parameter
nF
. Values
were chosen as acceptable fits if they fell within 0.3 log units of the global log10(
nF
) minimum.
Two examples of data-model fits – one closer to equilibrium, and one farther from equilibrium
– are shown in Figure S6. A contour plot of log
nF
versus
F
diss
and
r
f b
is shown in Figure
S6b. Panel
a
shows a profile of log
nF
as a function of
r
f b
at the best-fit value of log
R
diss
(red
line), and a profile of log
nF
as a function of
r
f b
along the covarying path of the local minimum
through model space (black line). Panel
d
shows
nF
across the entire model space. Panel
e
shows a profile of log
nF
as a function of
R
diss
at the best-fit value of
r
f b
(red line), and a profile
of log
nF
as a function of
R
diss
along the covarying path of the local minimum through model
space (black line). The best model fits typically had errors on the order of 10% or less, which are
comparable to the
∼
5% errors in experimental dissolution rate determinations (Subhas et al.
(2015)).
R
diss
and
r
f b
values were then binned based on the goodness of fit, and statistics were
collected on all
R
diss
and
r
f
b
pairs that were within 0.3 log units of the global
nF
minimum.
These best fits are found along the black-dashed line in Figure S6b. The distributions of
F
diss
and
r
f b
look very different for the two experiments in Figure S6. The high dissolution rate
experiment has a much longer tail in both
F
diss
and in
r
f b
, which indicates a large range in both
of these parameters that can adequately fit the data. Comparatively, the low undersaturation
experiment demonstrates a tighter distribution of
F
diss
and
r
f b
values, although there are a few
fliers at high dissolution rate that are paired with very low r(fb) values. We acknowledge that
the fit to the lower undersaturation experiment is systematically off, and there is still curvature
in this dataset that is not explained by the model. The suggestion of an underconstrained
saturation state is reasonable, but our experiments are designed to minimize changes in Ω
during a dissolution experiment. Both
δ
13
C and alkalinity provide measured constraints on the
magnitude of Ω before, during, and after a dissolution experiments. For high undersaturation
7
experiments where dissolution rate is fast, we run our experiment for less time, and monitor
in
situ
alkalinity changes due to dissolution. These changes result in at most a 0.03 change in Ω
over the course of a dissolution experiment. Thus we discount large changes in Ω affecting our
data from the perspective of our experimental design and our measurements of dissolution rate
and saturation state. Furthermore, the a change in Ω would affect the high dissolution rate data
more than the low dissolution rate data, and thus cannot be the source of all curvature in our
data sets, as near-equilibrium data sets are more curved than far from equilibrium data sets.
There is one other explanation for small curvature farther from equilibrium. The amount of
curvature expressed in the model is a function of the size of the fluxes relative to the size of the
reactive layer. At higher dissolution rate, the fluxes quickly reset the reactive layer and thus
very little curvature is expressed. Therefore, for this high net rate, there is very little curvature
in the model values that can fit this data set. One explanation for for the change in curvature
near Ω = 0
.
7 is that when the mechanism switches, the reactive layer deepens. Onset of 2D
nucleation across the entire mineral surface could indeed lead to a deepening of the reaction front
into the calcite lattice, as now the entire calcite surface is activated for dissolution. However,
as discussed in the main text, activation of the entire surface for dissolution also means that
precipitation reactions will have less of an influence on the isotopic composition of the solid, as
any precipitation should be immediately removed through re-dissolution. Given these caveats,
we conclude that the data fits using our model are valid, and continue with the interpretation
of our net rate data presented in the main text. All model files are available upon request to
the corresponding author.
Sensitivity to Various Model Parameters
Apart from
R
diss
and
r
f b
, the main model sensitivity is in the number of monolayers of calcite
(0.5 nm thick) that are allowed to react with solution. This parameter determines the shape
and strength of the curvature, and also affects the ability of the model to fit the dissolution rate
data. Figure S7 shows the goodness of fit for several different reactive layer thicknesses (1, 3, 5,
and 7 monolayers) as quartile box plots. The red line is the median of the misfit. Edges of the
boxes are the 25th and 75th percentiles. The whiskers include the data extremes and outliers
are plotted as red crosses. 5 monolayers provides the best overall misfit to the dataset. The box
model fits are slightly skewed: they are not evenly distributed around the mean misfit value.
This is partially because there is a floor at 0% misfit. But also, this skewness suggests that there
are high-misfit outliers. These higher misfits are due to
δ
13
C curvature that cannot be explained
using our fixed saturation-state box model and thus lead to a larger misfit in the
δ
13
C versus
time data. In Subhas et al. (2015), we calculated mean Ω values and rates for these experiments;
here, our model does not adjust its dissolution rate as a function of saturation state. We thus
attribute these highest misfits as a consequence of small changes in Ω, through either DIC loss or
alkalinity generation at high rates of dissolution, leading to a changing dissolution rate over the
course of the experiment. The misfit percentage does not scale with Ω, suggesting that our data
is not biased by an inability to sufficiently fit data as a function of saturation state. Based on the
overall ability of the model to fit all of our dissolution data, we chose a reactive layer thickness of
5 monolayers in our model, which was further justified through the SIMS experiments described
and discussed above and in the main text. More than 7 monolayers starts to become inconsistent
8
with the tracer incorporation measured via SIMS, as discussed above.
Expressions for the Dissolution Rate of Calcite
The canonical derivation of a linear dissolution rate law for calcium carbonate
Eq. 1 in the main text is a canonical representation of dissolution kinetics, and is the basis for
framing dissolution in terms of undersaturation, or 1-Ω:
R
net
=
R
f
−
R
b
=
k
f
{
CaCO
3
}−
k
b
[Ca
2+
][CO
2
−
3
] =
k
f
−
k
b
[Ca
2+
][CO
2
−
3
]
,
(10)
where
R
net
is in units of moles per time,
k
f
is in units of moles per time, and
k
b
is in units
of length
6
mole
−
1
time
−
1
. Assuming that the activity of the carbonate solid is 1,
R
f
=
k
f
and
R
b
=
k
b
[Ca
2+
][CO
2
−
3
]. Substitution of
K
′
sp
=
[Ca
2+
][CO
2
−
3
]
sat
{
CaCO
3
}
=
k
f
k
b
into Eq.(10), and the
definition of Ω =
[Ca
2+
][CO
2
−
3
]
K
′
sp
, gives:
R
diss
=
k
r
·
K
′
sp
−
k
r
[Ca
2+
][CO
2
−
3
] =
k
r
·
K
′
sp
(1
−
Ω) =
k
diss
(1
−
Ω)
.
(11)
Here,
k
diss
is in units of moles per time. Dividing Eq.(11) by surface area gives the specific
dissolution rate in units of moles length
−
2
time
−
1
. In addition, this model suggests that the
ratio of dissolution and precipitation rates is linearly proportional to Ω:
Ω =
[Ca
2+
][CO
2
−
3
]
K
′
sp
=
[Ca
2+
][CO
2
−
3
]
·
k
r
k
f
=
R
b
R
f
.
(12)
It is important to emphasize that curvature in either the gross rate of dissolution or precipi-
tation – and thus in the net rate of dissolution – as a function of saturation state necessitates
nonlinearities in either or both of these terms. For instance, curvature in the gross dissolution
rate implies that either the activity of the solid is not unity at all saturation states, or that there
are other terms that must be included in the rate law. The inclusion of other terms, for instance
carbonic acid as a driver of dissolution, may be the way forward. A fully predictive rate law
as a function of chemical specitaion would thus include multiple terms, and would necessitate
a speciation model of the calcite surface as a function of dissolved species and the density of
calcium and carbonate surface ion sites in seawater.
A derivation of Eq. 3 in the main text
In the last section, we discuss a 2D-nucleation model of calcite dissolution similar to that pre-
sented by Dove et al. (2005). This model takes elements of classical growth and nucleation
theory and applies them to dissolution. The rate of dissolution initiated by etch pit formation
is defined as:
R
n
=
hv
2
/
3
J
1
/
3
,
(13)
where the normal dissolution rate
R
n
is a function of the step height
h
, the speed of a moving
step
v
, and the steady-state etch pit nucleation rate,
J
. The nonlinear dependence of
R
n
9
on
J
is discussed extensively elsewhere (Sangwal, 1987; Malkin et al., 1989). The form of this
dependence changes the absolute value of surface energies calculated in Table 1 in the main text;
however, it does not impact the trends in surface energy or the location of the rate transition
between defect-initiated and homogeneous etch pit nucleation shown in Figure 5 in the main
text. The step retreat velocity
v
depends on solution composition in a formulation similar to
Eq.(10) above:
v
=
ωβ
(
C
e
−
C
) =
ωβC
e
(1
−
Ω);
(14)
where
ω
is the molecular volume of calcite,
β
is the step kinetic coefficient cm/s, and
C
e
is
the equilibrium concentration of the species. In the case of calcite,
C
e
is equivalent to the
K
′
sp
of calcite.
J
, the frequency of nucleating an etch pit, is related to saturation state
σ
= lnΩ,
as nucleation is activated above some ∆G
critical
, defined by the interface free energy barrier to
nucleation
α
:
J
=
|
σ
|
1
/
2
n
s
ahC
e
β
exp
[
πα
2
ωh
(
k
B
T
)
2
∣
∣
∣
∣
1
σ
∣
∣
∣
∣
]
.
(15)
New terms here are the lattice spacing
a
and the density of nucleation sites
n
s
. Substituting
Equations 15 and 14 into Eq.(13) and rearranging for 1
/σ
, we recover Equation 3 in the main
text.
Carbonic Anhydrase and its Effect on Dissolution Kinetics
The dissolution rate method from Subhas et al. (2015) allows for sensitive rate determinations
in the absence of significant changes of solution chemistry or mineral surface area. For instance,
using a 100% labeled CaCO
3
solid, we achieve a
δ
13
C sensitivity of about 20
h
per 1
μ
eq/kg
alkalinity change. Saturation state was determined using DIC-Alk pairs measured on a Picarro
CRDS and a home-built titration system. Final errors on Ω, calculated using Alkalinity-DIC
pairs, range from
≤
0.01 to
∼
0.03 units. The
K
′
sp
for calcite was adjusted so that the most
saturated dissolution experiments are undersaturated, requiring a correction factor of about
1.03 to the value found in Morse et al. (1980). The alkalinity contribution of carbonic anhydrase
(
∼
30 equivalents/mol) to solutions was determined by a standard additions alkalinity experiment
in natural seawater. At the same saturation state (Ω=0.83), curves of
δ
13
C increase significantly
in slope as CA increases (Figure S8).
To test the effect of other proteinaceous material on the rate of calcite dissolution, exper-
iments were conducted in the presence of Bovine Serum Albumin (BSA) at concentrations of
0.002 (not shown) and 0.01 mg/mL (shown in Figure S9 versus uncatalyzed and catalyzed disso-
lution rates at [CA] = 0.01 mg/mL). These experiments show no significant change in the BSA
rate versus the uncatalyzed rate while dissolution experiments in the presence of CA are always
faster than those without CA. If anything, rates in the presence of BSA are slightly slower.
The behavior of dissolution rate versus ∆G is very similar to that described by Arvidson and
L ̈uttge (2010). In their experiments, surfaces with different amounts of etching were exposed
to solutions of similar undersaturation. These surfaces dissolved at significantly different rates,
10
suggesting that the surface itself exerts significant influence on the relationship between disso-
lution rate and saturation state. Here, in contrast, it is the presence of CA, instead of a surface
feature hysteresis, driving a different functional dependence between ∆G and dissolution rate.
Figures and Tables
Profile Type
12
C
/
13
C
Mg/Ca
Beam type
Cesium
Oxygen
Sample HV
-5 kV
-8.5 kV
Aperture
300
μ
m
300
μ
m
Raster size
50x50
μ
m
50x50
μ
m
Beam current
0.5 nA
3 nA
Incidence Angle
24.5
◦
22.5
◦
Mass Resolution ∆m/m
3,000
2,000
Table S1:
Pertinent information for SIMS profile analysis on the CAMECA 7f instrument.
Ω = 0
.
95
Ω = 1
.
3
12
C
tot
0.7-2.2
·
10
−
13
2.1-6.3
·
10
−
13
Mg
2+
tot
1.5-7.8
·
10
−
15
1.3-8.8
·
10
−
15
z
ef f,
12
C
(nm)
1.0-3.3
3.1-9.3
z
ef f,Mg
(nm)
0.3-1.5
0.2-1.7
Table S2:
Estimates of new tracer incorporation based on integration of
12
C
/
13
C and Mg/Ca SIMS
profiles. The number of moles added were determined using Eq.(1) of SIMS profile data. This was then
converted to a thickness of calcite as described in the text.
11