Supplement of A model for marine sedimentary carbonate diagenesis and paleoclimate proxy signal tracking: IMP v1.0

Abstract. The preservation of calcium carbonate in marine sediments is central to controlling the alkalinity balance of the ocean and, hence, the ocean–atmosphere partitioning of CO2. To successfully address carbon cycle–climate dynamics on geologic (≫1 kyr) timescales, Earth system models then require an appropriate representation of the primary controls on CaCO3 preservation. At the same time, marine sedimentary carbonates represent a major archive of Earth history, as they have the potential to preserve how seawater chemistry, isotopic composition, and even properties of planktic and benthic ecosystems, change with time. However, changes in preservation and even chemical erosion of previously deposited CaCO3, along with the biogenic reworking of upper portions of sediments, whereby sediment particles are translocated both locally and nonlocally between different depths in the sediments, all act to distort the recorded signal. Numerical models can aid in recovering what the “true” environmental changes might have been, but only if they appropriately account for these processes. Building on a classical 1-D reaction-transport framework, we present a new diagenetic model – IMP (Implicit model of Multiple Particles (and diagenesis)) – that simulates biogeochemical transformations in carbonate-hosted proxy signals by allowing for populations of solid carbonate particles to possess different physicochemical characteristics such as isotopic value, solubility and particle size. The model also utilizes a variable transition matrix to implement different styles of bioturbation. We illustrate the utility of the model for deciphering past environmental changes using several hypothesized transitions of seawater proxies obscured by sediment mixing and chemical erosion. To facilitate the use of IMP, we provide the model in Fortran, MATLAB and Python versions. We described IMP with integration into Earth system models in mind, and we present the description of this coupling of IMP with the “cGENIE.muffin” model in a subsequent paper.


S1. Model options S1.1. Signal tracking method S1.1.1. Time-stepping method IMP can adopt the time-stepping method (method 1 in Fig. 3a) as an option. In this section, the time-stepping method is used for dissolution experiment #2 of Section 3.2.2 and results are 5 compared with those by method 2. The isotope signal changes in δ 13 C and δ 18 O are discretized into 42 steps, represented by 42 classes of CaCO 3 particles with distinct signals (Fig. S1) as described in Section 2.4.1.
When using the time-stepping method with 42 different CaCO 3 classes for dissolution experiment #2, very similar results are obtained to those with the interpolating method (com-10 pare Figs. S2a-c and S2d-f). Results of the time-stepping method approach the interpolation method results when the total number of CaCO 3 classes is increased. Furthermore, other simulations adopting the time-stepping method lead to the same degree of agreement of the results (e.g., dissolution experiment #1). Such comparisons confirm that signal tracking diagenesis can be simulated by either method 1 or 2 but a large number of CaCO 3 classes is required when 15 adopting method 1.

S1.1.2. Direct tracking of CaCO 3 isotopologues
As another option, IMP is able to directly track temporal changes of spatial distributions of five This option is a specific application of the direct tracking method (method 3 in Fig. 3c) and is beneficial because it allows tracking of 14 C age and ∆ 47 , in addition to δ 13 C and δ 18 O, in a computationally efficient way (Section 2.4.1). For 14 C age tracking the radioactive decay of Ca 14 CO 3 and accompanying generation of alkalinity (through assumed production of nitrate 14 N) are accounted for: where is one of n cc CaCO 3 classes including Ca 14 CO 3 and Ca 12 CO 3 , and λ is the decay constant which is 8033 −1 yr −1 for Ca 14 CO 3 (Aloisi et al., 2004) and zero for all other classes.
With method 3, Eqs. (S1) and (S2) can be implemented instead of Eqs. (5) and (9) for 14 C age 20 tracking. Apart from the implementation of Eqs. (S1) and (S2) instead of Eqs. (5) and (9) the model calculation and application to signal tracking are the same as described in Section 2.
As an example for the direct tracking of CaCO 3 isotopologues, we conduct again dissolution experiment #2 as presented in Section 3.2.2. Here changes in ∆ 47 in CaCO 3 are imposed additionally and are assumed synchronous with the δ 18 O changes (Fig. S3). The rain fluxes of 25 S1 individual isotopologues (Fig. S3b) are calculated directly from the input signals (Fig. S3a), except for the Ca 14 CO 3 flux, which is calculated assuming that the 14 C/ 12 C ratio is constant at 1.2 × 10 −12 (Aloisi et al., 2004). Note that the changes in rain fluxes with time are hard to discern from the semi-log plot of Fig. S3b. Another set of the five isotopologues needs to be additionally simulated when one further tracks model time as a proxy.

30
Obtained signals of δ 13 C and δ 18 O as well as CaCO 3 abundance (Fig. S4a-c) are the same as those obtained by method 2 (Figs. 13g-i), thus validating the direct tracking method. As the imposed temporal changes for ∆ 47 are the same as those for δ 18 O (Fig. S3a), the same patterns of signal distortion are obtained for both signals (Figs. S4b and e). The 14 C age is increased from ∼5 up to ∼50 kyr by dissolution (Fig. S4d), which is consistent with previous 35 results (e.g., Keir, 1984;Broecker et al., 1991;Oxburgh, 1998;Barker et al., 2007), but with our model, we can further examine the effects of different styles of bioturbation as well as different extents of dissolution (e.g., Fig. S4d). Note, however, that 14 C ages plotted against diagnosed depth do not reflect 14 C decay below the mixed layer (Section 2.4.2 and Fig. 5).
As an interesting feature of the option, we can further examine the kinetic isotope effect 40 on signal distortion. As an example, assuming −0.005% of kinetic isotope effect for dissolution To illustrate the utility of this model option, here we conduct a set of diagenetic experiments.
In the first, steady-state bulk profiles are obtained in the same way as we obtained the profiles in Fig. 6, except now the water depth is 4 km and the total sediment depth (z tot ) is 200 m ('CTRL' in Fig. S5). In the second run, we add 12 and 24 µmol cm −2 yr −1 of DIC and ALK 60 fluxes, respectively, at 10 m sediment depth, from assumed deep occurence of anoxic oxidation of methane (AOM) ('+AOM at 10 m' in Fig. S5). The two fluxes are added assuming normal S2 distributions centered at 10 m with 1 m standard deviation. We also allow precipitation of CaCO 3 by removing the Heaviside function from Eq. (5) for both experments.
Because of the addition of DIC and ALK fluxes from AOM, porewater gets supersaturated with respect to CaCO 3 at depth, and resultant authegenic CaCO 3 adds ∼20 wt% to the background CaCO 3 in the control run where AOM is not imposed. The two experiments demonstrate the utility of the option to examine the effect of a specific OM-related reaction on CaCO 3 diagenesis, aided also by the flexibility of the grid structure in IMP (compare Figs. S5 and 6; see also Section S1.4 CaCO 3 diagenesis and proxy signals.

S1.4. Grid structure of model sediment
The grid structure of modeled sediment is determined by the total sediment depth (z tot ) and a parameter β that determines how the sediment layer thickness changes with sediment depth: the grid becomes more irregular as β approaches 1 (Table 2; cf. Hoffman and Chiang, 2000).

85
The default grid structure adopted in the main text is a relatively deep and highly irregular sediment grid where z tot = 500 cm and β = 1 + 5 × 10 −11 . As a beneficial feature of IMP, one can adopt a different grid structure that suits the user's interest (cf. Section S1.2), although the changes in the above two parameters have to be made within the codes.
As an example of such a user option, we repeat the dissolution experiments in Section 3.2.2 90 but now adopting a shallower (z tot = 50 cm) and less irregular (β = 1.05) sediment grid. Despite significant difference in the grid structure, temporal changes of depth profiles are similar to those with the grid adopted in the main text (compare e.g., Figs. S8 and 7). Signal tracking methods in Section 2.4.2 can also be found efficient to minimize the effects of numerical diffusion on tracked signals regardless of the grid structure because very similar proxy signals are simulated S1.5. Time tracking We illustrated signal tracking diagenesis with the interpolating method (method 2, Section 2.4.1) in Section 3.2. Simulations in Sections 3.2.1 and 3.2.2 track 2 paleoceanographic proxies δ 13 C and δ 18 O with 4 classes of CaCO 3 particles that have the maximum and/or minimum values of 100 input δ 13 C and δ 18 O. In Section 3.2.3, we doubled the number of CaCO 3 classes (i.e., adopted 8 classes) to track the δ 13 C and δ 18 O in two distinctive model species ('fine' vs. 'coarse' species).
By adopting 8 classes of CaCO 3 particles that have the endmember values of δ 13 C and δ 18 O and size-dependent dissolution and mixing properties, one can track not only δ 13 C and δ 18 O, but also the size distribution in bulk CaCO 3 in simulations in Section 3.2.3 (cf. Fig. 15c). Thus, a 105 simple relationship can apply to method 2: when one tracks n p proxies and/or characteristics As examples, the following subsections (Sections S1.5.1-S1.5.3) repeat the same simulations 120 as in Sections 3.2.1-3.2.3 but with enabling tracking of model time by adopting 8 or 16 classes of CaCO 3 particles whose properties are tabulated in Tables S1 and S2. Except for the doubled numbers of simulated CaCO 3 classes (Tables S1 and S2), the experimental conditions are the same as those presented in Section 3.2.

125
The experiments to examine distortions of proxy signals by bioturbation in Section 3.2.1 are repeated in this section with additionally tracking model time (  Table S1. Before describing the results, the details on the calculation of the rain fluxes for individual CaCO 3 classes are presented first. When one considers a given proxy or characteristic, CaCO 3 130 classes can be grouped into 2 sets of classes, which have either the maximum or minimum value for the cosidered proxy/characteristic (Table S1). The rain fractions for the two sets can be calculated from the input value for the proxy/characteristic (e.g., Fig. S10b when considering δ 13 C). For instance, if the two sets of CaCO 3 classes have the values of A and B and the input value is given as X (cf. Fig. 3b), then the rain fractions for the two sets of classes (denoted 135 as a and b, respectively) can be calculated as a solution of two equations: aA + bB = X and a + b = 1. After repeating the assignment of the rain fractions for two sets of classes for all the considered proxies/characteristics (Figs. S10a-d), the rain fraction of a given class can be calculated as the product of the rain fractions for the sets that have the same properties as those of the considered class (Fig. S10e). As an example, rain fraction of class #1 (dotted On the other hand, a single age model cannot be applied to reveal temporal changes of 150 proxy signals when carbonate dissolution is mechanistically linked to styles of bio-mixing (Figs. S11g and h). Apparent signal shifts in depth caused by changes in the extent of dissolution or sediment accumulation rate (Figs. S11e and f) are absent when plotted against mode time using individual age models (Figs. S12d and e). It also becomes clearer that proxy signals are more effectively lost along with CaCO 3 particles when mixing is relatively ineffective. For example, a 155 slight reduction of δ 18 O peak by enabling dissolution is recognized for the no bioturbation case (see a slightly smaller peak of pink curve in Fig. S12e than in Fig. S12b).

S1.5.2. Dissolution of carbonates
The same dissolution experiments as in Section 3.2.2 ( Fig. 12) are repeated in this section except with enabling model-time tracking. Instead of 4 we adopt 8 classes (Table S1) to track 160 model time in addition to δ 13 C and δ 18 O. The temporal changes of the rain fractions for the 8 CaCO 3 classes are presented in Fig. S13b, calculated in a way described in Section S1.5.1 (Table S1 and Fig. S13b).
Perfect age models for different cases with and without bioturbation are available with enabling time tracking, with which signal change events recorded in the sediment (i.e., Fig.   165 13) can be plotted against model time (Fig. S14) S15).
Using the tracked time against diagnosed depth (i.e., age model) proxy signals in the two distinctive CaCO 3 species can be plotted against model time (Fig. S16). As described in Section 3.2.3, apparent offsets of excursions of proxy signals and species abundance between the two 180 species disappear (compare Figs. S16 and 15) but the reduction of the magnitude of imposed species abundance changes is still recognized especially in cases with bioturbation (Fig. S16c).          In a-c, dissolution rate constants of all CaCO 3 classes are fixed at zero, while in d-f, at the default value (Table 1). S19 Figure