Seasonal cycles of the carbon export flux in the ocean: insights from the SISSOMA mechanistic model

This study aims to investigate the seasonal cycle of carbon export flux in the ocean using the SISSOMA modeling framework. SISSOMA uses a 2-dimensional state space (size and excess density) to follow the fate of aggregates in the mixed layer which are transformed through three main processes, e.g., aggregation, fragmentation, and remineralization, until they eventually sink out of the surface ocean. The model tracks aggregate size, mass, and porosity which allows for a direct estimate of aggregate sinking speed through a Reynolds number modified Stokes' law. First, a simple seasonal cycle with a single peak of POM production is presented, which provides a solid basis to understand the model's framework and enables us to perform sensitivity analyses on important parameters. The effect of increased stratification on a reconstructed ecosystem in the north Atlantic is then presented and discussed. Overall, our results showcase the nonlinear relationship between the production of primary particles and the export of aggregates out of the mixed layer and unveil key mechanics of the three transformative processes. Moreover, it has been shown that remineralization rates, stickiness, and the size/ excess density characteristics of the primary particles all affect in various ways the intensity, seasonal cycle, and the resulting size spectrum of the aggregate community. Finally, our results indicate the crucial role that turbulence plays in both the timing and the magnitude of the carbon export flux which might affect not only the potential of the system to remove carbon out of the mixed layer but also have a direct impact on the organisms inhabiting the mesopelagic layer which rely on the sinking particles to cover their energetic needs.

1 Introduction

The production and export of particulate organic matter (POM) from the surface ocean to depth constitutes an important pathway in the biogeochemical cycle of the ocean, maintaining in part the carbon stored in the deep ocean by the biological carbon pump (BCP) (Boyd et al., 2019; DeVries, 2022), as well as supplying benthic and mesopelagic ecosystems with organic carbon and energy (Billett et al., 1983; Hernández-León et al., 2020). It is qualitatively well understood that export is ultimately governed by primary productivity in the sunlit surface ocean. Although a mechanistic description of how these are related remains elusive (Boyd and Newton, 1995; Buesseler, 1998; Henson et al., 2012), a simple heuristic measure, the e-ratio (the ratio of export production F to net primary production P_NPP, sometimes referred to as export efficiency) is often invoked in models and analyses (Siegel et al., 2014; Henson et al., 2019; Buesseler et al., 2020; Laurenceau-Cornec et al., 2023). Factors that contribute to the variation in e-ratio are myriad and include the structure of the plankton community (Boyd and Newton, 1995; Guidi et al., 2016; Henson et al., 2019), nutrient cycling (Kemp et al., 2000; Raven and Waite, 2004), turbulence and cell exudates (Passow et al., 2001; Engel, 2004).

A somewhat simplifying perspective can be achieved by dividing the e-ratio into two components (Laufkötter et al., 2016); the first being the relationship between P_NPP and P_POM which is the production rate of primary particulate organic material . This component (P_POM $/$ P_NPP), termed the p-ratio (Laufkötter et al., 2016), is largely controlled by the plankton community and is generally the most easily accessible component of the export flux simulated in ecosystem models. The second component is the relationship between P_POM and export, F. This component is largely controlled by aggregation processes, transforming primary detrital material into fast sinking aggregates, and it is here that perhaps the greatest uncertainty lies in estimating export flux. The overall estimation of the e-ratio can thus be broken down as:

$$\begin{array}{}\text{(1)}& e\text{-ratio}={\displaystyle \frac{F}{P_{\mathrm{NPP}}}}={\displaystyle \frac{P_{\mathrm{POM}}}{P_{\mathrm{NPP}}}}{\displaystyle \frac{F}{P_{\mathrm{POM}}}}=p\text{-ratio}\times s\text{-ratio}\end{array}$$

(Laufkötter et al., 2016). It is the second factor in this relationship (termed the s-ratio; Laufkötter et al., 2016) that is the focus of this study. Notably, recent models seeking to resolve the e-ratio (Siegel et al., 2014; Nowicki et al., 2022; Laurenceau-Cornec et al., 2023) from satellite products do not consider aggregation per se, but rely on empirically derived sinking speeds for different classes of material i.e. focusing almost exclusively on the p-ratio. While doubtless there are important ecological drivers contained in the p-ratio, the aggregation processes in the s-ratio could exert an equally strong control from time to time.

Here we explore the seasonal cycle of POM production and export using the SISSOMA modeling framework. Given the production rate and characteristics (e.g. size, ballasting, mass) of primary POM particles, we simulate the formation and transformation of POM into aggregates while estimating sinking speed directly via Stokes' law corrected for a finite Reynolds number. We simulate not only aggregation, but also remineralization and fragmentation, two important but often poorly resolved aspects of aggregate dynamics. In particular, the sinking velocity of the resulting aggregates is crucial for their fate, as fast-sinking aggregates are subjected to lower remineralization losses by bacterial in the surface and eventually reach greater depths in the oceans' interior. The aggregate's sinking velocity depends on its size, excess density, composition (e.g. ratio of organic matter to ballast minerals) and porosity, while remineralization rates are mainly defined by the aggregate's lability (food quality), the structure of the microbial community and temperature (Baumas and Bizic, 2024). Finally, it is well understood that the relationship between P_POM and F exhibits a strong nonlinearity (Kiørboe et al., 1994). This is qualitatively described as a critical concentration (Jackson, 1990; Burd and Jackson, 2009) where aggregation to larger fast sinking aggregates accelerates rapidly and can lead to episodic intense export events (Siegel et al., 2025) precipitated by intensifying turbulence. It is far from clear that any simple relationship exists between P_POM and F, and by extension F and P_NPP. This is reflected in the high variability in the observed e-ratio (Henson et al., 2012, 2019; Siegel et al., 2023) both spatially and temporally, suggesting that a more mechanistic analysis is required.

2 Methodology

2.1 Model description

SISSOMA is a mechanistic model (developed in MATLAB) that tracks the fate of aggregates through three main processes; aggregation, fragmentation, and remineralization, and their subsequent sinking. These processes are simulated in a 2-dimensional state space specifying size and excess density of the aggregates. The model accounts for variations of aggregate porosity, mass, and size, allowing for an estimate of the aggregate sinking speed while they are constantly being redistributed in the 2-dimensional state-space through transformative processes. As input, the aggregation model can use estimates of the production/input of primary particles, representing for instance dead and dying plankton cells, fecal pellet production, and/or aolean dust deposits, each characterized by their size and excess density.

Many of the key characteristics of primary particles can be linked to their origin, for instance diatoms have a higher excess density than other phytoplankton taxa (Hansen and Visser, 2019), and the size of the fecal pellets are governed by the size of zooplankton grazers (Stamieszkin et al., 2015). For the purposes of this study, we construct primary particle production rates that can be inferred from typical seasonal cycles of the plankton community. All the descriptions and values of the variables used for the current analysis can be found in Table E1 in Appendix E.

A key feature of SISSOMA is that it treats aggregation as a self-similar process and introduces a novel convolution algorithm for application in a multidimensional state space. This is in contrast to other aggregate models that either use a Monte Carlo approach (Jokulsdottir and Archer, 2016) and/or rely on simplifying aggregate descriptors such as uniform fractal dimension or particle size spectrum slope. Importantly, fractal dimension cannot be considered a conservative property under remineralization and fragmentation, a shortcoming in existing models that SISSOMA seeks to address.

The physical space representation of SISSOMA is kept as simple as possible: a uniform surface mixed layer throughout the annual cycle, although this can be expanded to multiple layers stacked to each other. The state space represents the size (radius) r and excess density ρ^′, (ρ−ρ_w), of the emerging aggregate community. SISSOMA provides flexibility to choose the range and the number of classes of this space.

It is convenient to introduce two transformed variables (x,y) that map to $(r,\mathit{\rho }-{\mathit{\rho }}_{\mathrm{w}})$ as:

$$\begin{array}{}\text{(2)}& r=r_o{\mathit{\delta }}^x,\mathit{\rho }-{\mathit{\rho }}_{\mathrm{w}}={\mathit{\rho }}_o{\mathit{\lambda }}^y{\mathit{\delta }}^{(a-\mathrm{3})x}\end{array}$$

where a is the self-similarity parameter, x is a logarithmic scaling factor of aggregate size, y a stretched logarithmic scaling factor of excess density, $\mathit{\delta }=(r_{\max }/r_o{)}^{\mathrm{1}/X}$ is the logarithmic radius interval, and $\mathit{\lambda }=({\mathit{\rho }}_{\max }/{\mathit{\rho }}_o{)}^{\mathrm{1}/Y}$ is the excess density logarithmic interval. The factor δ^(a−3)x takes advantage of the reduction of density by aggregation and expands the density resolution for large aggregates. r_o, ρ_o and r_max, ρ_max are the minimum and maximum values of the aggregate's size and density, respectively.

For this project, the size of possible aggregates in the system, r, ranges between 1 and 10⁶ µm and consists of 30 logarithmically spaced bins, while their excess density, ρ^′, extends between 1 and 100 kg m⁻³, Fig. 1a. Throughout this domain, sinking speed w, Fig. 1b, is estimated via a modified version of Stokes' law:

$$\begin{array}{}\text{(3)}& w^{\mathrm{2}}={\displaystyle \frac{\mathrm{8}}{\mathrm{3}}}{\displaystyle \frac{{\mathit{\rho }}^{\prime }}{{\mathit{\rho }}_{\mathrm{w}}}}{\displaystyle \frac{gr}{C}}\end{array}$$

where ρ_w is the density of the surrounding water, C is the drag coefficient (see Table E1), and g is the gravitational acceleration.

Figure 1 Representation of the (a) aggregate excess density, ρ^′, and (b) sinking speed, w, in the 2-dimensional state space where x is a logarithmic scaling factor of aggregate size, y a stretched logarithmic scaling factor of excess density (see Eq. 2). (a) Larger aggregates lead to lower excess densities, as surrounding water is trapped in their interior during the formation period. The colors of the dashed boxes in (a) point out the three scenarios in the sensitivity analysis in Sect. 3.2, where different size/excess density primary particles are introduced in the system: (i) small and light (S + L) in grey, small and dense (S + D) in green, and big and dense (B + D) in purple. (b) The sinking speed of aggregates as derived from Eq. (3).

A rigorous mathematical description of particle aggregation dates back more than a century (Smoluchowski, 1918). Following more recent formulations (Burd and Jackson, 2009; Jokulsdottir and Archer, 2016) and including specific representations of fragmentation and degradation, aggregate dynamics can be summarized as follows.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}N(\mathit{s},t)& =\\ & q_{m:c}{P}_{\mathrm{POM}}^{\prime }(\mathit{s},t)\\ & +{\displaystyle \frac{\mathit{\alpha }}{\mathrm{2}}}\underset{\mathrm{S}}{\int }\mathit{\beta }({\mathit{s}}^{\prime },c^{\star }(\mathit{s},{\mathit{s}}^{\prime }\left)\right)N\left(c^{\star }\right(\mathit{s},{\mathit{s}}^{\prime }),t)N({\mathit{s}}^{\prime },t)\mathrm{d}{\mathit{s}}^{\prime }\\ & -\mathit{\alpha }N(\mathit{s},t)\underset{\mathrm{S}}{\int }\mathit{\beta }(\mathit{s},{\mathit{s}}^{\prime })N({\mathit{s}}^{\prime },t)\mathrm{d}{\mathit{s}}^{\prime }\\ & -\mathit{\psi }\left(\mathit{s}\right)N(\mathit{s},t)+\underset{\mathrm{S}}{\int }\mathit{\psi }\left({\mathit{s}}^{\prime }\right)N(s^{\prime },t)\mathit{\theta }(s,s^{\prime })\mathrm{d}{s}^{\prime }\\ & -{\displaystyle \frac{\partial \mathit{s}}{\partial t}}\cdot {\displaystyle \frac{\partial N}{\partial \mathit{s}}}-N(\mathit{s},t){\displaystyle \frac{w\left(\mathit{s}\right)}{h}}\end{array}\end{array}$$

where N(s) is the number of aggregates in a given state-space bin, $\mathit{s}=(r,{\mathit{\rho }}^{\prime })$ is a state variable which includes relevant aggregate attributes and the function c^⋆ defines the state of a parent aggregate based on the state of both the other parent and the produced-daughter aggregate. Primary particles are produced in the mixed layer at a rate of $q_{m:c}{P}_{\mathrm{POM}}^{\prime }$, where q_m : c represents the dry mass to carbon ratio, and the second term is the time and state-resolved production rate of POM. That is, $P_{\mathrm{POM}}^{\prime }\mathrm{d}\mathit{s}$ is the production rate of POM dry mass in the same interval of state space $[\mathit{s},\mathit{s}+\mathrm{d}\mathit{s}]$. Note that $P_{\mathrm{POM}}={\int }_{\mathrm{S}}{P}_{\mathrm{POM}}^{\prime }\mathrm{d}\mathit{s}$, i.e. the integral over all state space.

The aggregation process is captured by the following two terms and relies both on the encounter kernel of two particles, $\mathit{\beta }(\mathit{s},{\mathit{s}}^{\prime })$, and their stickiness, α, which defines the probability of two particles to form an aggregate once they have collided, see Appendix B. The larger an aggregate, the more prone it is to break down at a rate of ψ(s), and to be redistributed to smaller aggregates with variable excess density according to the partitioning function $\mathit{\theta }(\mathit{s},{\mathit{s}}^{\prime })$, see Appendix D. The next term describes loses from remineralization, which we assume occurs at a constant rate throughout the state space and removes only dry mass leading to aggregates of the same size-class but reduced excess density, see Appendix C. Finally, whether aggregates are exported out of the mixed layer, h, depends on their sinking speed, w(s).

Key variables in the model are the matrices N and M representing the number of aggregates and their total dry mass respectively within bins in discretized (x,z) state space. Suitable range choices for $x\in [\mathrm{0},X]$ and $y\in [\mathrm{0},Y]$, scaling factors r_o [µm] and ρ_o [kg m⁻³], and logarithmic intervals δ and λ allow for a relatively complete representation of the aggregate community within computationally convenient dimensions of N and M. Total aggregate dry mass M and number density N are related by $M=m^{\prime }\circ N$ where ∘ represents piece-wise matrix product and m^′ is the mean dry mass of an aggregate within each bin.

2.2 Analyses description

In the first part of the analysis, a highly seasonal environment is reproduced with a cycle of primary particle production peaking in the middle of April at 1000 mgC m⁻² d⁻¹ , Fig. 2a. To better understand the model output, given a set of parameters (see Table E1), we distribute the primary particles along the smallest size class of 1 µm and distribute them equally throughout the excess density field, Sect. 3.1. Next, we investigate the effect of remineralization, stickiness and primary particles' size/excess density characteristics on the export flux, Sect. 3.2. For this purpose, we assign the primary particles along a size range between 1 to 50 µm and distribute them equally throughout the excess density field, Fig. 2b. Regarding the latter, we simulate three scenarios based on the relative characteristics of the primary particles: (a) small and light, S + L (in grey), (b) small and dense, S + D (in green), and (c) big and dense, B + D (in purple) primary particles, as shown in Fig. 1a. Further, we investigate the effect of episodic storm events on export flux with regards to their intensity (i.e. the level of turbulent dissipation rate), their duration, and their occurrence in different times of the seasonal cycle, see Sect. 3.2.

Figure 2 Reconstruction of (a) a simple annual production cycle of primary particles peaking in middle April, (b) which are assigned to a size range between 1–50 µm and are equally distributed throughout the excess density space. x is the logarithmic scaling factor of aggregate size, y is the stretched logarithmic scaling factor of excess density (see Eq. 2). This scenario is the base for the results in Sect. 3.2, 3.3 and 3.1 with the exception that in the latter analysis only primary particles of size 1 µm are produced in this case study.

In the second part of the analysis, a unicellular plankton community in the north Atlantic is replicated where a strong “spring bloom” and a weaker “autumn bloom” are observed. In this case, we distinguish between primary particles that originate from diatoms and the rest of the unicellular organisms (i.e. bacteria, phytoplankton, heterotrophic microzooplankton and mixotrophs). For simplicity, we term the latter as generalists (Cadier et al., 2020; Serra-Pompei et al., 2020), Fig. 3a. In the beginning of the “spring bloom”, the optimal light, temperature and stratification conditions, as well as the high availability of nutrients, give a competitive advantage to diatoms. They grow faster until a point where nutrients (especially silicate) become limited and the rest of the unicellular organisms take over as their smaller sizes enables them to take up nutrients more efficiently at low concentrations. The same pattern but in lesser magnitude is reproduced in the case of the “autumn bloom”. It is important to mention that, in the context of this project, in the seasonally stratified scenario only the turbulent dissipation rates varies with time, while the mixed layer depth remains constant at 100 m throughout the annual cycle in all simulations. Regarding the distribution of the primary particles in the 2-dimensional space, we assume that when an organism dies it produces POM of the same size. Diatoms' radius ranges between 3.2 and 50.6 µm, whereas the rest of the unicellular organisms' radius varies between 1 and 32 µm. Moreover, based on the fact that diatoms posses a silica shell and their excess density is higher than the rest of the unicellular community, we distribute the primary particles originating from diatoms in the upper half of the excess density space while primary particles introduced by the rest of the unicellular organisms are equally allocated throughout the full length of the excess density space, Fig. 3b. Finally, we investigate how the mixing conditions in the upper 100 m of this environment affect the transformative processes of POM (namely aggregation and fragmentation) and the resulting carbon export. This is accomplished by focusing on the turbulent dissipation rate through two scenarios: one with high seasonal variability and a second scenario with a relatively constant rate, Fig. 3c. In both scenarios, the mixed layer depth is not affected by these changes and stays constant at 100 m in time.

Figure 3 Reconstruction of an idealized ecosystem in the north Atlantic with primary materials originating from two functional groups, i.e., diatoms and generalists, used for the analysis in Sect. 3.3. (a) The annual cycle is characterized by two distinct “bloom” periods where new material enter the system, (b) primary particles originating from diatoms are assigned to larger size bins and are equally distributed to higher excess density ranges compared to generalists (x is the logarithmic scaling factor of aggregate size, y is the stretched logarithmic scaling factor of excess density, see Eq. 2), (c) two scenarios were used to represent the stratification conditions of the system: a highly seasonally-mixed and a stratified throughout the year. Both scenarios are based on a well-mixed, vertically uniform water column layer with a constant depth of 100 m, where the turbulent dissipation rate is constant in depth but variable in time.

Throughout this analysis, three different annotations are used to describe the export flux (mgC m⁻² d⁻¹): f₁₀₀ refers to the density integrated export flux (Eq. 5), F₁₀₀ is the size and excess density integrated export flux in each point in time (Eq. 6) and F_tot is the total, annually integrated export flux.

$$\begin{array}{}\text{(5)}& {\displaystyle }{f}_{\mathrm{100}}\left(x\right)=\sum _y{m}_{\mathrm{dry}}(x,y)w(x,y)\text{(6)}& {\displaystyle }{F}_{\mathrm{100}}=\sum _{x,y}{m}_{\mathrm{dry}}(x,y)w(x,y)\end{array}$$

where m_dry is the total dry mass concentration (mgC m⁻³) and w is the sinking velocity (m d⁻¹). Finally, following (Laufkötter et al., 2016) we define instantaneous s-ratio as the fraction of the size and excess density integrated export of material in each time (F₁₀₀) to the production of new, primary particles (P_POM) in that time:

$$\begin{array}{}\text{(7)}& s\text{-ratio}={\displaystyle \frac{F_{\mathrm{100}}}{P_{\mathrm{POM}}}}\end{array}$$

3 Results

3.1 Model analysis

In the first part of the results, the analysis on a simple seasonal cycle with primary particles of radius 1 µm produced in the upper mixed layer (see Sect. 2.2) is presented. Figure 4a shows how both the production of new material (P_POM) and the export flux (F₁₀₀) evolve in time, as well as illustrates the relationship between them expressed as the instantaneous s-ratio. A time-lag of 68 d between the peaks of P_POM and F₁₀₀ is observed, while the s-ratio fluctuates between 0.001 and 0.65.

Figure 4 A simple annual cycle where the smallest size-class (1 µm) primary particles with variable excess densities are produced in the mixed layer. (a) The progression of the production (P_POM), total export flux (F₁₀₀) and instantaneous s-ratio in time (blue solid line), (b) the distribution of dry mass in the size-spectrum in each time-step, (c) the annual evolution of the density-integrated export flux, (f₁₀₀), and (d) the system-characteristic phase diagram where the total export flux (F₁₀₀) in time is plotted against the POM production (P_POM). The black dot in (d) highlights the annual mean export flux and production resulting in an annual mean s-ratio of 0.13.

During the first four months of the annual cycle, there is an accumulation of primary particles and production of small proto-aggregates in the system, Fig. 4b. This does not yet lead to a strong export event which can be explained by the fact that these first aggregates are small, in the range of 1 to 200 µm, and sink very slowly. However, there is a distinct shift in the size-spectrum of the dry mass in the system where aggregates with almost five orders of magnitude of difference in their size range are created, Fig. 4b. This implies that there is a critical concentration of mass at which the aggregation process becomes the main transformative force, producing aggregates of variable, bigger sizes and excess density which in turn are able to support the following strong export event, Fig. 4c, peaking at a total of 182 mgC m⁻² d⁻¹.

Finally, the phase diagram in Fig. 4d highlights the strong non-linear relationship between the production of new material and the resulting export flux. The black dot indicates the annual mean fluxes of production and export which results in an annual mean s-ratio of 0.13.

3.2 Sensitivity analysis

In the second part of the analysis of this system, a set of sensitivity analysis on three parameters; remineralization, stickiness and size/excess density characteristics of primary particles, is presented, Fig. 5. Figure 5a, b show the effect of remineralization on the export flux and its emergent size-spectrum, respectively. The higher the remineralization rate, the less material stays in the system over time to later being transformed through aggregation, fragmentation, remineralization. This is directly reflected on the magnitude of the resulting export flux. Higher remineralization rates result in weaker export flux as reflected in the phase-diagram, the annual mean s-ratios in Table 1 and the annually integrated export fluxes in the enclosed graph in Fig. 5b.

Figure 5 Phase diagram (a, c, e) and density-integrated export flux (f₁₀₀) (b, d, f) for a range of different (a, b) remineralization rates γ (day⁻¹), (c, d) stickiness α (–) and (e, f) size excess density characteristics of the primary particles (S + L refers to small and light, S + D to small and dense and B + D to big and dense). The dots denote the annual mean export flux and production for each scenario. Regarding the density integrated export flux figures, the shadings represent the annual variability, whereas the lines with the circle marks are the annual mean values. The enclosed figures in (b), (d), (g) show the annually integrated export flux (mgC m⁻² yr⁻¹) for each scenario.

Table 1 Annual mean s-ratios for different remineralization and stickiness values and primary particles characterstics

The stickiness of aggregates affects export flux in a reversed way compared to remineralization, Fig. 5c, d. The annual mean s-ratio increases with increasing stickiness, Table 1, which can be related to the effectiveness of the system to export carbon out of the mixed layer. In all cases, the dominant size classes of the resulting size spectrum are around 10⁴ µm with the difference that higher stickiness results in higher variability in export during the annual cycle, Fig. 5d.

The size/excess density characteristics of the primary particles also play a crucial role in the development of the export flux over time, Fig. 5e, f. Our analysis shows that big and dense primary particles (B + D) lead to the highest total export flux annually, resulting in an annual mean s-ratio of 0.87, Table 1. Moreover, the phase-diagram indicates that once these B + D primary are produced, they are either exported immediately or very quickly form bigger (fast-sinking) aggregates. However, in the case of small and dense (S + D) primary particles, there is a time-lag of approximately 60 d between the start of the production cycle and the increase in export flux. Comparing these two scenarios, a system with B + D primary reaches its peak export flux earlier in the year (day 125 compared to day 171 for the S + D scenario), and its size-spectrum is dominated by smaller aggregates, Fig. 5f. Finally, small and light (S + L) primary particles cannot support significant export flux, as the majority of the production is lost due to remineralization before reaching the critical mass concentration and form faster-sinking aggregates.

Figure 6 Representation of the effect of the size and excess density of the primary particles on the export efficiency of the system. F₁₀₀ is the total export flux and P_POM is the production of new (primary) particles in the system. The loop (phase diagram) shows the annual cycle of the system. The red dot denotes the start of cycle (1 January) and the loop progresses anti-clockwise. The three simulated sizes (radius) are 1.27, 31.9 and 801.4 µm. Given these sizes, the respective examined excess densities are [1.6 10.2 64], [0.06 0.4 2.6] and [0.002 0.02 0.1] kg m⁻³. The slope of the dotted lines (see values) indicate the annual mean ratio of F₁₀₀ to P_POM (s-ratio).

Figure 6 presents the effect of the characteristics of the primary particles, e.g. size and excess density, in a systematic way. As one would expect, a combination of larger and heavier primary particles leads to the highest export efficiency among all scenarios. In this case, primary particles can either sink fast in their one or form even faster-sinking aggregates in a short time period, which in both cases minimizes the losses due to remineralization. In the other extreme, small and light primary particles need a lot of time to form larger, fast-sinking aggregates. During this period, remineralization is the dominant transformative process, removing material from the system in a fast rate and hence minimizing its export efficiency. An interesting insight of this analysis is that, given a specific size of primary particles, increasing their excess density level has a significant and positive impact on the system's export efficiency. However, given a specific excess density level of primary particles, the effect of increasing size is suggested to plateau faster. This figure also showcases how responsive the export flux is to the primary particle production, which is reflected both in the shape of the phase diagram and the range of the instantaneous ratio of export flux to the production of primary particles (see the colorbar). For example, focusing on the case of the highest excess density of 63 kg m⁻³ but the smallest primary particles, there is a great mismatch between production and export with the highest s-ratio values recorded at the end of the production cycle.

3.3 The effect of turbulence

Turbulence is another important variable as it is involved in both the aggregation and fragmentation processes. Higher turbulence, on the one hand, increases the encounter rate of two particles and, on the other hand, it makes large aggregates more vulnerable to break down into smaller aggregates of variable excess density. Sudden and strong storm events have the potential to mix the water column and increase its turbulence levels. Figure 7 attempts to simulate such storms varying both their duration and intensity. Figure 7a shows that the longer the duration of the storm (at the peak of the primary particle production), the stronger the export flux, with the effect being more pronounced for higher turbulent dissipation rates in each scenario. When comparing these scenarios, we can also identify a distinct shift around August. Before that point, the longer storm events produce the stronger export flux whereas the opposite pattern is observed (even though it is less pronounced) after this point. This might be an indication that these strong episodic events, by removing large amounts of material out of the mixed layer, inhibit the system maintaining the critical mass concentration needed for aggregation in the following months. On the other hand, a system with a more conservative and steady behavior keeps its material in the system for longer time, where even though a part of it is lost through remineralization, the remaining material keeps being repackaged into bigger aggregates. This can be better visualized in Fig. 7b where during a 14 d storm event in early September, the scenario of 1 d storm event in spring responses more strongly to increased turbulence. An interesting behavior of the system is observed during the second storm event in Fig. 7b, where at the beginning of the storm the export flux decreases before rising in a steep slope to reach its peak. This implies that fragmentation is at first the dominant process, filling the system with smaller-slow sinking aggregates, but providing the necessary “fuel” for the aggregation processes later on.

Figure 7 (a) Simulation of sudden storm events of variable duration (t_storm) and intensity (reflected by the turbulent dissipation rate, ϵ). The blue-dotted line shows the POM production cycle (P_POM) whereas the colored lines demonstrate the export flux (F₁₀₀) in each scenario. The enclosed graphs show the annually integrated export flux (mgC m⁻² yr⁻¹) for each scenario. The background turbulent dissipation rate is $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{6}}$ m² s⁻³ and the turbulent dissipation rate of the “fall” storm event in (b) is $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{4}}$ m² s⁻³ and lasts for 14 d in all scenarios.

Figure 8 The effect of turbulence on the export flux, as analyzed in a highly seasonal (a, c) and a rather constant turbulent environment (b, d). The idealized ecosystem reconstruction of this analysis includes two types of unicellular organisms producing primary particles of different size/excess density characteristics. The mixed layer depth is constant at 100 m in both cases and only the turbulent dissipation rate varies in time as indicated in Sect. 2.2 and Fig. 3a, b present the size-spectrum of export flux (density integrated flux) and (c, d) show the phase diagram with the annual mean export flux and production being denoted in the black dot for each scenario respectively.

In the final part of the analysis, the effect of turbulence in a highly seasonal environment with primary particles originating from an imitated unicellular plankton community is investigated, see Sect. 3.3. SISSOMA's output suggests that higher magnitude and variability of turbulence throughout the year leads to higher total export fluxes, i.e., 3.02×10⁴ mgC m⁻² yr⁻¹ compared to 1.92×10⁴ mgC m⁻² yr⁻¹ in the scenario of a rather steady and low turbulent dissipation rates, Fig. 8. Hence, the efficiency of the system to export the produced particles out of the mixed layer varies with an annual mean s-ratio of 0.44 in the former case and 0.28 in the latter. Higher turbulence might imply that the system possesses a lower critical mass concentration point (where aggregation boosts the creation of bigger, fast-sinking aggregates) which in turn can be reached faster. This can be observed by comparing the evolution of the export flux in time, Fig. 8a, b, where the highly seasonal system supports a stronger flux earlier in the year and hence in lower concentrations of material in the system. On the other hand, the delay observed in the low turbulence scenario leaves primary particles and slow-sinking speed aggregates suspended in the water column and exposed to remineralization which decreases the potential of the system to export that material out of the mixed layer. This can be better visualized in the phase diagram, Fig. 8c, d, where during the autumn bloom the export flux in the low turbulence scenario collapses, even though the input of new material in the system is the same for both scenarios.

4 Discussion

SISSOMA provides a modeling framework to mechanistically describe the formation and export of aggregates out of the mixed layer. It approaches this by incorporating information about both aggregate size and excess density which are constantly transformed in this 2-dimensional state space, through three main processes, e.g., aggregation, fragmentation, remineralization, leading to their eventual sinking. The current work focused on how these processes play out over a seasonal cycle and investigated how a variety of parameters, such as remineralization rate, stickiness, size/excess density characteristics of the primary particles and turbulence, all affect the intensity and the size structure of the resulting export flux over time. The results of sensitivity analyses underscore the need for more observational and experimental studies on how microbial degradation acts on marine particles and finding a systematic way to represent how the stickiness of different sourced marine particles varies during the seasonal cycle. In addition, the importance of gaining information about the size/density characteristics of the primary particles entering the mixed layer has been uncovered. Regarding the computational cost of SISSOMA, in its current form, it is more suitable to be applied for understanding the export flux on a regional scale or testing its results against experimental setups. Future work with the help of machine learning techniques may help us to generalize SISSOMA's architecture, which could enable its use in global biogeochemical models. This requires the deployment of SISSOMA simulations under a wide range of primary particle characteristics (size/ excess density) and parameters (such as stickiness), the relation of these results to different ecosystem types in the global ocean, and finally the acquisition of enough in situ data for training and testing the performance of the machine learning algorithm. In this direction, a generalized and computationally-light POM specialized model will enable us to more accurately describe the fate of POM in the upper mixed layer, nutrient availability dynamics and thus plankton community structure when integrated in more complex global biogeochemical models.

First, our results highlight the strong nonlinear relationship between the production and the export of material out of the mixed layer, and how these can be described and quantified via the mechanisms involved in the transformative processes. Aggregation is sensitive to the concentration of mass in the system, as our analysis has shown there is a clear transition where the formation rate of bigger aggregates accelerates rapidy, a behavior that supports the concept of “critical cell concentration” (Jackson, 1990). These larger, fast-sinking aggregates not only are exported rapidly but also avoid the high remineralization losses experienced by slower sinking aggregates. This threshold of accelerated aggregation, in turn, is dependent on the size/excess density characteristics of the primary particles, their stickiness and the prevailing turbulence in the system. In the first part of the analysis, the observed increase in export flux is very sharp and sudden which might appear strange, but it can be explained by the fact that the new-primary particles produced in the mixed layer are assigned to the smallest size-class (1 µm). In general, the rate of aggregation depends on three encounter kernels, each governed by relative aggregate motion induced by Brownian motion, turbulence and the differential settling, respectively. By assuming a constant background turbulence, the aggregation potential of small particles depends on the Brownian motion which means that the aggregation process is very slow, as the smaller the particles the lower the chances encountering each other (Burd and Jackson, 2009). During this time, material accumulates in the mixed layer at a rate determined by production and remineralization until a point where the concentration is high enough for aggregation to dominate, a process that is very sudden and accelerates rapidly. The story is different when primary particles of variable (and higher) size-classes enter the system, Sect. 3.2. In this case, both differential settling and turbulent shear contribute to a faster and more efficient encounter of primary particles or proto-aggregates which leads to a more direct response of the export flux to the introduction of new material in the mixed layer. This is also reflected on the instantaneous and annual mean s-ratios, which are tracking the different export efficiencies between the two systems. The large range of fluctuations of the instantaneous s-ratio throughout the year suggests that care should be taken when generalizing the results of short-term field expeditions. The measured efficiency of carbon export strongly depends on the sampling time and delay effects between production and export of material. In this direction, we could even use SISSOMA's modeled annual s-ratio cycle to apply corrections on the measurements and draw more representative conclusions. Overall, the size/excess density characteristics of the primary particles directly affect the mass-concentration threshold for the aggregation to become the dominant transformative process which in turn defines the observed time-lag between production and export flux. As POM is an important energy source for a wide array of organisms, especially in the mesopelagic zone, even small changes in timing might lead to a mismatch between their phenology and the availability of their food with unpredictable cascading effects (Robinson et al., 2010).

Moreover, the smaller the primary particles, the further they are from the size/ excess density ranges which can lead to a fast sinking out of the mixed layer and every time a new aggregate is formed their relative excessive density decreases (due to the inclusion of surrounding water). In a system like this, this implies that a portion of aggregates might end up becoming big but not dense enough to sink. Once these aggregates reach the size-threshold where fragmentation is dominant, the 2-dimensional state-space fills with aggregates of various sizes and excess densities. In the other extreme scenario of big and dense primary particles entering the system, this effect of fragmentation in time is not so strong. These primary particles can form aggregates which are either rapidly exported or their relatively high density fragments are repackaged and get exported in short time. Understanding these mechanisms is very important when we want to predict how shifts in the structure of the plankton community (and hence shifts in the characteristics of primary particles) will affect the export fluxes in the future, with direct consequences on the global carbon cycle and the availability of energy sources in the deeper ocean.

Our results indicate that the level of this “critical mass concentration” is also dependent on turbulence and stickiness in a way that higher values of these variables lead to a lower critical mass concentration point. This implies that even in the smallest mass concentrations, particles have higher encounter rates and stick to each other more efficiently moving mass into bigger, fast-sinking size classes. The first part of the analysis in Sect. 3.3 shows how short, episodic events of high turbulence, such as during wind storms, affect the evolution of the export flux in time and its magnitude. A strong storm event might at that specific moment enhance the export of material out of the system, however, this leads to a reduced capability of the system to export material in the following period, as their concentration in the mixed layer stays low and aggregation processes weaken. Moreover, our results indicate that the duration of these events defines the magnitude of their impact on the export flux. When comparing a highly seasonally mixed system to a stratified system (even though this was a comparison of two extreme situations with the turbulent dissipation rate varying, but the export depth being kept constant), our results motivate us to look closer into the consequences that increased stratification in deep convection regions, such as the north Atlantic, could have in their potential to export and sequester carbon. Overall, the seasonally mixed system supports a stronger export flux which comes earlier in the year, whereas the stratified scenario failed to sustain an “autumn” export event. Apart from all the above-mentioned processes, fragmentation plays its role in a way that depending on which point of the annual mixing cycle the aggregate fragments are released in the system, they can either rapidly be repackaged and exported out of the mixed layer (in the case of high turbulence) or remineralized in case the turbulence stays low for long period.

Finally, seeing holistically the export system with all three transformative processes in play, we could identify three characteristic stages, Fig. 9. First, there is a “preparation phase” where new primary particles are produced in the mixed layer, a portion of which might sink fast out of it (based on their size/excess density properties), but the rest of the material accumulates in the system. By remineralization being the dominant process, a portion of them is lost to nutrients and progressively moves into lower excess density bins. This ends when a critical mass concentration level is reached where aggregation becomes the dominant process, moving mass to larger size bins and in intermediate excess density bins compared to the “parent” particles or proto-aggregates. The stickiness, turbulence and the relative proportion of size/density primary particle characteristics all affect the behavior of the system and the critical level in which each phase takes place. It is more difficult to draw certain conclusions about the transition between the export and de-escalation phases. However, according to our understanding higher stickiness, turbulence and proportion of big/dense primary particles might lead to a shorter (but more intense) export phase draining the system out of necessary mass for aggregation to keep being a dominant transformative process. Fragmentation constantly redistributes mass to smaller size classes of various excess densities. The third and last phase starts when the mass concentration falls below the critical level. Now, aggregation and fragmentation continue to operate, but remineralization's influence increases steadily. In this way, all three transformative processes are in play, but the system is progressively occupied by lower size and importantly lower excess density material. These slower-sinking aggregates stay now suspended in the water column and most of them are lost due to remineralization, until new production comes and the cycle starts again. The timing of these events and how fast they occur in this cycle depends on both the size/excess density characteristics of the primary and the environmental conditions, as discussed above.

Figure 9 Summary of an annual cycle of the export system as explored in SISSOMA. The width of the arrows and the radius of the circles represent the relative importance of the three transformative processes. The arrows point towards the direction where mass is forced to move by each process. Three distinct period can be identified: (1) a preparation phase where remineralization is the dominant transformative process acting on the primary particles and the small proto-aggregates, (2) the main “export phase” is established whenever the critical mass concentration, C_crit, is reached. This depends on turbulence, stickiness and the size/ excess density characteristics of the primary particles. The inset plot shows their effect on the C_crit in an abstract way, where the x-axis refers to a spectrum of possible values of each parameter of interest. During this period, aggregation, fragmentation and sinking out of the mixed layer are the dominant processes redistributing mass in the 2-dimensional state space and eventual exporting a big part of it, (3) when the system has exported a big portion of the material in the mixed layer (the critical mass concentration falls below a threshold) the “de-escalation” phase starts where all processes are in play until the cycle starts again.

Even though SISSOMA captures many important factors and processes related to an ecosystem's potential to export material out of the mixed layer, there are certain assumptions that need to be taken into account and improved. In its current form, remineralization removes dry mass from an aggregate in a constant rate throughout the year and moves it to lower excess density bins. SISSOMA can provide a platform for testing a variety of hypotheses on how biological, physical and chemical factors, as well as the aggregates' properties all affect the degradation rates of the latter. Anderson et al. (2023) suggest that microbial remineralization roughens the surface area of the sinking aggregates which in turn increases the attachment area for the microbes to colonize and eventually degrade it faster than expected. In our case, the incorporation of porosity as a third state variable for the aggregate's characteristics could allow us to incorporate this information when more experimental data are set to draw a mathematical formulation of this process. Moreover, modeling the effect of ballast minerals on sinking velocities and degradation rates is crucial. Dense mineral particles, such as calcium carbonate or opal, usually become an important part of an aggregate increasing its sinking velocity. This, on the one hand, shortens the period that this newly-formed aggregate stays in the upper ocean layers where remineralization is higher. On the other hand, there is contradicting evidence that ballast minerals can create a protective layer around an aggregate, further decreasing the microbial degradation rates (Cram et al., 2018). It has also been suggested that disproportionally high ballast mineral concentrations might have totally opposite overall effect by producing more but smaller aggregates, due to the low proportion of the necessary organic matter – “glue” for keeping them together in the long run (De La Rocha et al., 2008). Finally, the way particles are moved around the 2-dimensional state space during their degradation time can be challenged, too. Contrary to the current direction of mass, one might argue that by removing the organic-light part of an aggregate, the proportion of ballast-heavy minerals inside the aggregate increases which in this case might result in heavier, fast-sinking aggregates instead.

Regarding the prevailing environmental conditions of the system, we assumed a uniformly mixed layer of constant depth throughout the year. However, the situation in nature is more variable (especially in high latitudes) with wide fluctuations in the strength of stratification and the depth of the evolved mixed layer, which might have an immediate effect on the export flux (de Boyer Montégut et al., 2004). A shallower mixed layer depth after a phytoplankton bloom, for example, might restrict the produced material into a confined area which could increase their collision rate and hence transform them into big, fast-sinking aggregates. However, rapid changes in stratification might move material out of the mixed layer, a process known as the mixed layer pump (Dall'Olmo et al., 2016), reducing the amount of material available in the mixed layer which then affects the relative strength between transformative processes.

Another relevant property to test in the model is the variable lability of the primary particles. It has been proposed that new production initiated by the input of excess nutrients, mainly in upwelling regions and deep convection areas, is labile and an easy target for the microbes which turn a big proportion of the sinking material into nutrients. In contrast, highly recycled material in the tropics and oligotrophic regions is usually refractory (Francois et al., 2002). This might imply that small, suspended material or bigger but not heavy enough aggregates remain in the system for longer periods, form aggregates of higher sinking speed and support considerable export fluxes in time. In the current SISSOMA version, the stickiness of the primary particles remains constant throughout the annual cycle. Grønning and Kiørboe (2022) showed that when certain diatom species sense the presence of copepods in their close proximity, they can quickly increase their stickiness to form big, fast sinking colonies to escape predation. Moreover, it has been suggested that during a phytoplankton bloom particles of varying stickiness are produced depending on a wide array of biological, physical and chemical conditions of the ecosystem (Baumas and Bizic, 2024). Even though we acknowledge the importance of stickiness which is also supported by our results, there is a need of further research in the topic to be able to confidently draw the appropriate mathematical framework and incorporate it in models. The above-described research areas could be approached by programming separate classes for primary particles coming from different sources, which then the model follows throughout the simulation. In this way, we could be able to know the proportions of different material in an aggregate, e.g., ballast minerals, fecal pellets, phytoplankton cells, in any time and decide more efficiently the direction of mass movement.

Last but not least, it has been suggested that the structure and phenology of the plankton community are the main factors that determine the characteristics of the export flux (Henson et al., 2019). In this direction, the next step is to couple SISSOMA with the Nutrient-Unicellular-Multicellular (NUM) model, a trait-based framework which uses first principles to model the global biogeography of the plankton community (Andersen and Visser, 2023).

5 Conclusion

In conclusion, SISSOMA provides a useful tool for mechanistically describing the seasonal cycle of carbon export flux. This framework helps us to better understand how aggregation, fragmentation and remineralization shape the emergent aggregate community in time, until their eventual export out of the mixed layer. Our results highlight the nonlinear relationship between the production of primary particles in the upper mixed layer and their export as aggregates out of it, which can be reflected on the wide fluctuations of instantaneous s-ratio throughout the annual cycle. Moreover, it has been shown how remineralization rates, stickiness, size/excess density characteristics of the primary particles and turbulence all affect both the intensity and the size-structure of the export flux. Although the theory and implementation of the aggregation processes are well established, more work needs to be done to improve the function of the fragmentation and remineralization processes in the model. This stresses the need for more, both in situ and experimental, research to be conducted, which will help us develop the appropriate mathematical framework and incorporate them into more sophisticated biogeochemical models.

Appendix A: Self-similarity

SISSOMA assumes that aggregation is a geometrically self-similar process, such that the linear dimension r_i∘j of an aggregate formed by the combination of two parent aggregates of linear dimension r_i and r_j respectively is given by:

$$\begin{array}{}\text{(A1)}& r_{i\circ j}=(r_i^{\mathrm{a}}+r_j^{\mathrm{a}}{)}^{(\mathrm{1}/a)}\end{array}$$

That is, for the binary process of aggregation, r^a is an additive conservative property. This is not a new idea (Jackson, 1998; Wiesner, 1992) and arises from the general observation that aggregates display fractal characteristics (Alldredge and Gotschalk, 1988; Meakin, 1987; Logan and Wilkinson, 1990). We term a the self-similarity parameter, and note a<3 in compliance with the observed increase in porosity under aggregation. In a departure from other aggregation descriptions, we do not interpret a as the fractal dimension of aggregates. Neither is it an inherent property of aggregates and we will not attempt to use a to produce scaling laws as is usual in these considerations (e.g. DeVries et al., 2014; Kriest and Evans, 1999; Cael et al., 2021). Specifcially while the self-similarity of aggregation is conserved, under the action of remineralization and fragmentation, the fractal dimension of aggregates is not a conservative property.

Under geometric self-similarity, the mass of an aggregate produced by the combination of 2 aggregates of mass m_i and m_j can be deduced to be the sum of these two masses, plus a bit extra due to the inclusion of some fluid (density ρ_w) that occupies the expanded aggregate volume (increased porosity). Specifically,

$$\begin{array}{}\text{(A2)}& m_{i\circ j}=m_i+m_j+{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }(r_{i\circ j}^{\mathrm{3}}-r_i^{\mathrm{3}}-r_j^{\mathrm{3}}){\mathit{\rho }}_{\mathrm{w}}\end{array}$$

Here m is the total mass, not just the dry mass of the aggregate. It is convenient to recast this in terms of density of the aggregates, ρ_i, ρ_j, ρ_i∘j, and their respective volumes, v_i, v_j and v_i∘j. It follows that excess density

$$\begin{array}{}\text{(A3)}& {\mathit{\rho }}_{i\circ j}-{\mathit{\rho }}_{\mathrm{w}}={\displaystyle \frac{v_i}{v_{i\circ j}}}({\mathit{\rho }}_i-{\mathit{\rho }}_{\mathrm{w}})+{\displaystyle \frac{v_j}{v_{i\circ j}}}({\mathit{\rho }}_j-{\mathit{\rho }}_{\mathrm{w}})\end{array}$$

That is, aggregation produces a daughter aggregate whose excess density is the volume weighted mean of the excess density of the two parent aggregates. Equations (A1) and (A3) provide an iterative computational framework that can follow the size and excess density and hence sinking speed, of aggregates as they are formed and transformed from primary particles.

Appendix B: Aggregation algorithm

Aggregation is computationally the most complex aspect of the model as it involves an off-set binary convolution of N in state space. Performing binary convolution calculations is greatly facilitated by self-similarity. Specifically, the ordinates of an aggregate produced from the combination of (x_i,y_i) and (x_j,y_j) is given by:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{x}_{i\circ j}& ={\displaystyle \frac{\mathrm{1}}{a\ln \mathit{\delta }}}\ln \left({\mathit{\delta }}^{a{x}_i}+{\mathit{\delta }}^{a{x}_j}\right)\\ y_{i\circ j}& ={\displaystyle \frac{\mathrm{1}}{\ln \mathit{\lambda }}}\ln \left({\mathit{\lambda }}^{y_i}{\mathit{\delta }}^{a(x_i-x_{i\circ j})}+{\mathit{\lambda }}^{y_j}{\mathit{\delta }}^{a(x_j-x_{i\circ j})}\right)\end{array}\end{array}$$

It follows that the combination of aggregates from any two domains of identical shape in (x,y) state space will be confined to a third identically shaped domain. The model utilizes this feature on a regular (x,y) grid in optimizing the algorithm architecture. Without loss of generality, aggregate state space can be defined in terms of transformed coordinates as

$$\begin{array}{}\text{(B2)}& \mathit{s}=\left[\begin{array}{c}x\\ y\end{array}\right]\end{array}$$

from which it follows that the binary combination of aggregates can be written as

$$\begin{array}{}\text{(B3)}& {\mathit{s}}_{i\circ j}=c^{\circ }(s_i,s_j)=c^{\circ }\left(\left[\begin{array}{c}{x}_i\\ y_{\mathrm{1}}\end{array}\right],\left[\begin{array}{c}{x}_j\\ y_j\end{array}\right]\right)=\left[\begin{array}{c}{x}_{i\circ j}\\ y_{i\circ j}\end{array}\right]\end{array}$$

From Eq. (B1), it can be easily shown that

$$\begin{array}{}\text{(B4)}& s_{i\circ j}^{\prime }=c^{\circ }({\mathit{s}}_i+s^{\prime },{\mathit{s}}_j+s^{\prime })={\mathit{s}}_{i\circ j}+s^{\prime }\end{array}$$

where s^′ is an offset in state space. This means that any bounded region in state space convolved with another bounded region of the same shape will map to a third bounded region again of the same shape. Notably, within discretized state space, aggregation from two 1×1 bins maps to a third 1×1 bin (Fig. B1).

Figure B1 State-space aggregation. Aggregates formed by interactions between bins s_i and s_j will be found in the region s_i∘j which will span, at most 4 bins in discretized state space.

For the binary aggregation ${\mathit{s}}_{i\circ j}\leftarrow c^{\circ }({\mathit{s}}_i,{\mathit{s}}_j)$, the time rate of change of particle numbers in appropriate state space bins are given by

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}\mathrm{d}N\left({\mathit{s}}_i\right)& =-\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_i\right)\\ \mathrm{d}N\left({\mathit{s}}_j\right)& =-\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_j\right)\\ \mathrm{d}N\left({\mathit{s}}_{i\circ j}^{\mathrm{00}}\right)& =\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_i\right)N\left({\mathit{s}}_j\right){\mathit{\xi }}^{\mathrm{00}}\\ \mathrm{d}N\left({\mathit{s}}_{i\circ j}^{\mathrm{01}}\right)& =\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_i\right)N\left({\mathit{s}}_j\right){\mathit{\xi }}^{\mathrm{01}}\\ \mathrm{d}N\left({\mathit{s}}_{i\circ j}^{\mathrm{10}}\right)& =\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_i\right)N\left({\mathit{s}}_j\right){\mathit{\xi }}^{\mathrm{10}}\\ \mathrm{d}N\left({\mathit{s}}_{i\circ j}^{\mathrm{11}}\right)& =\mathit{\alpha }\mathit{\beta }({\mathit{s}}_i,{\mathit{s}}_j)N\left({\mathit{s}}_i\right)N\left({\mathit{s}}_j\right){\mathit{\xi }}^{\mathrm{11}}\end{array}\end{array}$$

where

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\mathit{s}}_{i\circ j}^{\mathrm{00}}& =\left[\begin{array}{c}{x}_{i\circ j}^{\mathrm{00}}\\ y_{i\circ j}^{\mathrm{00}}\end{array}\right]=\left[\begin{array}{c}\mathrm{floor}\left(x_{i\circ j}\right)\\ \mathrm{floor}\left(y_{i\circ j}\right)\end{array}\right]\\ {\mathit{s}}_{i\circ j}^{\mathrm{10}}& ={\mathit{s}}_{i\circ j}^{\mathrm{00}}+\left[\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right]\\ {\mathit{s}}_{i\circ j}^{\mathrm{01}}& ={\mathit{s}}_{i\circ j}^{\mathrm{00}}+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right]\\ {\mathit{s}}_{i\circ j}^{\mathrm{11}}& ={\mathit{s}}_{i\circ j}^{\mathrm{00}}+\left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\end{array}\right]\end{array}\end{array}$$

and

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{\mathit{\xi }}^{\mathrm{00}}& =(x_{i\circ j}-x_{i\circ j}^{\mathrm{00}})(y_{i\circ j}-y_{i\circ j}^{\mathrm{00}})\\ {\mathit{\xi }}^{\mathrm{10}}& =(\mathrm{1}-x_{i\circ j}+x_{i\circ j}^{\mathrm{00}})(y_{i\circ j}-y_{i\circ j}^{\mathrm{00}})\\ {\mathit{\xi }}^{\mathrm{01}}& =(x_{i\circ j}-x_{i\circ j}^{\mathrm{00}})(\mathrm{1}-y_{i\circ j}+y_{i\circ j}^{\mathrm{00}})\\ {\mathit{\xi }}^{\mathrm{11}}& =(\mathrm{1}-x_{i\circ j}+x_{i\circ j}^{\mathrm{00}})(\mathrm{1}-y_{i\circ j}+y_{i\circ j}^{\mathrm{00}})\end{array}\end{array}$$

The total number of states s in discretized state space is L=XY, from which it follows that the number of unique binary combinations is $K=(L+\mathrm{1})L/\mathrm{2}$. For a moderate 10×8 state space resolution, K=3240. Defining $k=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }K-\mathrm{1}$ as an ordinate counting combinations, $\mathrm{\ell }=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }L-\mathrm{1}$ as an ordinate counting discrete states s, and

$$\begin{array}{}\text{(B8)}& \mathit{\kappa }\left(k\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{2}L+\mathrm{1}-\sqrt{(\mathrm{2}L+\mathrm{1}{)}^{\mathrm{2}}-\mathrm{8}k}\right)\end{array}$$

it follows that

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{\mathrm{\ell }}_i& =\mathrm{floor}\left(\mathit{\kappa }\right(k\left)\right)\\ {\mathrm{\ell }}_j& =k-{\mathrm{\ell }}_{i}L+{\mathrm{\ell }}_i({\mathrm{\ell }}_i+\mathrm{1})/\mathrm{2}\end{array}\end{array}$$

Further, each ℓ can be decomposed into x,y ordinates as:

$$\begin{array}{}\text{(B10)}& x=\mathrm{floor}(\mathrm{\ell }/Z);y=\mathrm{\ell }-xY\end{array}$$

That is, each k is associated with a unique combination $({\mathrm{\ell }}_i,{\mathrm{\ell }}_j)\leftrightarrow (x_i,y_i,x_j,y_j)$ from which the offset ordinates $({\mathrm{\ell }}_{i\circ j}^{\mathrm{00}},{\mathrm{\ell }}_{i\circ j}^{\mathrm{01}},$ ${\mathrm{\ell }}_{i\circ j}^{\mathrm{10}},{\mathrm{\ell }}_{i\circ j}^{\mathrm{11}})$ for binary aggregation can be defined as in Eq. (B6). These provide a means of mapping indices for binary combinations, to state space and back. A consequence of the model architecture is that once a state space is defined, all coefficients determining the convolution process (encounter kernels β, offset indices ${\mathrm{\ell }}_{i\circ j}^{\mathrm{00}},{\mathrm{\ell }}_{i\circ j}^{\mathrm{01}},{\mathrm{\ell }}_{i\circ j}^{\mathrm{10}},{\mathrm{\ell }}_{i\circ j}^{\mathrm{11}}$, and partitioning functions ${\mathit{\xi }}_{i\circ j}^{\mathrm{00}},{\mathit{\xi }}_{i\circ j}^{\mathrm{00}},{\mathit{\xi }}_{i\circ j}^{\mathrm{00}},{\mathit{\xi }}_{i\circ j}^{\mathrm{00}}$ ) are fixed. That is, these need only be calculated once, stored in vectors of length K, and referenced via a lookup map.

Appendix C: Remineralization

For any aggregate we can express its density ρ as

$$\begin{array}{}\text{(C1)}& \mathit{\rho }={\displaystyle \frac{m}{v}}={\displaystyle \frac{m_{\mathrm{d}}+m_{\mathrm{w}}}{v}}=(\mathrm{1}-p){\mathit{\rho }}_{\mathrm{d}}+p{\mathit{\rho }}_{\mathrm{w}}\end{array}$$

where p is the porosity, ρ_d is the density of the solid phase of the aggregate and ρ_w is the density of water. This leads to the expression for excess density ρ^′

$$\begin{array}{}\text{(C2)}& {\mathit{\rho }}^{\prime }=\mathit{\rho }-{\mathit{\rho }}_{\mathrm{w}}=\mathit{\varphi }({\mathit{\rho }}_{\mathrm{d}}-{\mathit{\rho }}_{\mathrm{w}})\end{array}$$

where we have used $\mathit{\varphi }=\mathrm{1}-p$. The simplifying assumption we make in this work is that within a short time step, reminealization removes dry mass from the aggregate, while the overall volume of the aggregate remains invariant. The volume of dry mass that is lost, is replaced by water thus increasing aggregate porosity p. The fundamental degradation process can be written as

$$\begin{array}{}\text{(C3)}& {\displaystyle \frac{\mathrm{d}{m}_{\mathrm{d}}}{\mathrm{d}t}}=-\mathit{\gamma }{m}_{\mathrm{d}}\end{array}$$

Writing m_d=vϕρ_d and assuming that the density of the dry mass phase remains invariant it follows that

$$\begin{array}{}\text{(C4)}& {\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}t}}v{\mathit{\rho }}_{\mathrm{d}}=-\mathit{\gamma }\mathit{\varphi }v{\mathit{\rho }}_{\mathrm{d}}\end{array}$$

or more succinctly

$$\begin{array}{}\text{(C5)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}^{\prime }}{\mathrm{d}t}}=-\mathit{\gamma }{\mathit{\rho }}^{\prime }\end{array}$$

Thus, from the expression for ρ^′ in terms of ordinates x,y (Eq. 2)

$$\begin{array}{}\text{(C6)}& {\displaystyle \frac{\mathrm{d}{\mathit{\rho }}^{\prime }}{\mathrm{d}t}}={\mathit{\rho }}^{\prime }\ln \left(\mathit{\lambda }\right){\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}=-\mathit{\gamma }{\mathit{\rho }}^{\prime }\end{array}$$

That is

$$\begin{array}{}\text{(C7)}& {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}=-{\displaystyle \frac{\mathit{\gamma }}{\ln \left(\mathit{\lambda }\right)}}\end{array}$$

Appendix D: Fragmentation

Fragmentation is simulated simply as a rate at which aggregate mass is transported to smaller sizes classes. We implement this as an increasing function of aggregate size. Further, it is intuitively understood that fragmentation increases as turbulent stresses increase. Of the processes considered, fragmentation remains the least well constrained; aggregates appear resistant to mechanical shear (Alldredge et al., 1990), and fragmentation appears to be chiefly mediated by metazooans through handling and feeding appears to be important (Dilling and Alldredge, 2000) and by microbial “mining” and dissolution of adhesive material. We set the fragmentation rate at ψ₀=0.5 d⁻¹ for large porous aggregates, a value consistent with observations (Briggs et al., 2020). Specifically, we set

$$\begin{array}{}\text{(D1)}& \mathit{\psi }(x,y)={\mathit{\psi }}_{\mathrm{0}}r(\mathrm{1}-\mathit{\varphi }){\displaystyle \frac{\mathit{\epsilon }}{{\mathit{\epsilon }}_c}}={\mathit{\psi }}_{\mathrm{0}}{r}_o{\mathit{\delta }}^x\left(\mathrm{1}-{\mathit{\delta }}^{(a-\mathrm{3})}\right){\displaystyle \frac{\mathit{\epsilon }}{{\mathit{\epsilon }}_c}}\end{array}$$

where ψ₀ is the fragmentation rate, r is the radius, (1−ϕ) is an estimation of the aggregate's porosity, $\mathit{\delta }=(r_{\max }/r_o{)}^{\mathrm{1}/X}$ is the radius logarithmic interval, a is the self-similarity parameter and ε_c is the reference turbulence. Given the uncertainty in fragmentation processes, we choose a simple redistribution rule ($\mathit{\theta }(\mathit{s},{\mathit{s}}^{\prime })$ in Eq. 4) where the mass from fragmented aggregates is simply equally distributed across all smaller size classes. Specifically, the mass of fragments $\mathrm{d}M(x,y)=\mathit{\psi }\left(x\right)M(x,y)$ from $\mathit{s}=(x,y)$ is distributed equally across all bins $\mathit{s}=(x^{\prime },y)$ for $\mathrm{0}>x^{\prime }>x$.

Appendix E

Figure E1 The three components of the coagulation kernels: (a) brownian motion, (b) differential settling, and (c) turbulent shear.

Table E1 Glossary of variables, parameter settings for simulations and units.