A Stochastic , Lagrangian Model of Sinking Biogenic Aggregates in the Ocean ( SLAMS 1 . 0 ) : Model Formulation , Validation and Sensitivity

We present a new mechanistic model, Stochastic Lagrangian Aggregate Model of Sinking particles (SLAMS) for the biological pump in the ocean, which tracks the evolution of individual particles as they aggregate, disaggregate, sink, and are altered by chemical and biological pro5 cesses. SLAMS considers the impacts of ballasting by mineral phases, binding of aggregates by transparent exopolymer particles (TEP), zooplankton grazing, and the fractal geometry (porosity) of the aggregates. Parameterizations for age-dependent organic carbon (orgC) degradation kinetics, 10 and disaggregation driven by zooplankton grazing and TEP degradation, are motivated by observed particle fluxes and size spectra throughout the water column. The model is able to explain observed variations in orgC export efficiency and rain ratio from the euphotic zone and to the sea floor as driven 15 by sea surface temperature and the primary production rate and seasonality of primary production. The model provides a new mechanistic framework with which to predict future changes on the flux attenuation of orgC in response to climate change forcing. 20


Introduction
Plankton in the ocean incorporate dissolved carbon and nutrients into particulate form (cells, marine snow and fecal pellets), which allows the constituents to sink through the water column until the particles decompose at depth.This process, called the biological pump (Ducklow et al., 2001;Buesseler et al., 2008), results in a major rearrangement of the chemistry of the oceans.The CO 2 concentration in the atmosphere is coupled with the chemistry of the surface ocean, giving the biological pump a hand in controlling atmospheric CO 2 and the climate of the Earth (Schneider et al., 2008;Kwon et al., 2009).One major uncertainty in future climate change projections is the response of the biological pump in the ocean to surface ocean conditions (thermal and chemical) (Friedlingstein et al., 2001).Here we develop a new numerical model, stochastic, Lagrangian aggregate model of sinking particles (SLAMS), to simulate the processes and characteristics of sinking particles in the ocean, to gain a better understanding of what factors affect the flux of organic carbon (orgC), and how it might respond to changing upper ocean chemistry and climate.
Particles in the ocean span about 5 orders of magnitude in size and 15 orders of magnitude in number density (Stemmann and Boss, 2012).The slope of the particle-size spectrum is highly variable (Guidi et al., 2009).Measurements have shown that the majority of the flux is composed of slow (< 10 m day −1 ) and fast (> 350 m day −1 ) sinking particles (Riley et al., 2012;Alonso-Gonzalez et al., 2010).
Flux attenuation of sinking organic particles is determined by the sinking velocity of particles and the rate of orgC respiration.Sinking velocity generally increases with particle size and density (Alldredge and Gotschalk, 1988;Ploug et al., 2008) and it has also been found to increase with depth, driven by evolution of particle size or density (Berelson, 2002).Normalized Particulate organic carbon (POC) flux to mass flux shows that orgC comprises about 5 % of the mass flux in the deep ocean (Armstrong et al., 2002).The near constancy of this ratio suggests that ballasting by minerals plays a key role in determining the orgC sinking flux (Klaas and Archer, 2002;Francois et al., 2002;Boyd and Trull, 2007).On the other hand, minerals are not inherently sticky, and an overabundance of mineral particles relative to transparent exopolymer particles (TEP) in laboratory experiments has been Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Jokulsdottir and D. Archer: Sinking aggregates found to decrease the sizes of the aggregates (Passow and De La Rocha, 2006) and sinking velocity (Passow et al., 2012), which could act to diminish the sinking flux.
Biologically and chemically driven processes act on particles to change their physical and chemical structure.The respiration of organic matter by bacteria is a complex, sequential process (Biddanda and Pomeroy, 1988), as indicated by an observed correlation between the degradation rate of orgC and its age, spanning 8 orders of magnitude (Emerson and Hedges, 1988;Hedges and Keil, 1995).Association with mineral surfaces also acts to protect orgC from enzymatic degradation (Engel et al., 2009;Le Moigne et al., 2013).Three mechanisms take place in zooplankton guts: respiration of orgC, dissolution of calcium carbonate and repackaging of phytoplankton cells into fecal pellets.Fecal pellets that are packaged more tightly than aggregated marine snow have a potential to sink quickly through the water column.Aggregate fragmentation is driven by zooplankton swimming and feeding (Alldredge et al., 1990) and the bacterial breakdown of TEP, which leaves aggregates un-sticky and prone to break up (Passow and De La Rocha, 2006).In addition to biological dissolution of CaCO 3 in zooplankton guts, there appears to be significant dissolution of the more soluble CaCO 3 mineral aragonite in the water column (Gangstøet al., 2008).Biogenic silica, the other main biomineral in the ocean, is undersaturated throughout the ocean, and is distinguished by a strong temperature sensitivity of its dissolution rate, dissolving significantly faster in warm waters (Gnanadesikan, 1999).
A common approach to modeling the flux of material through the water column in large-scale ocean models is using the Martin curve (Martin et al., 1987), which predicts the orgC flux to depths from the export, F (z 0 ), and an attenuation parameter, b: Its simple form makes it well suited for large models where computational economy is important, such as Ocean Carboncycle Model Intercomparison Project -phase II (OCMIP-II) (Najjar et al., 2007).The difficulty lies in understanding how b varies regionally and temporally (Guidi et al., 2009(Guidi et al., , 2015)).The model of Armstrong et al. (2002) accounts for the ballast effect, that minerals affect the orgC flux by protecting it from degradation or increasing its sinking velocity, by modeling two classes of orgC: one that is associated with ballast minerals and one that is not.Community Climate System Modelocean Biogeochemical Elemental Cycle (CCSM-BEC) expands on the ballast model to simulate particle fluxes and accounts for temperature and O 2 to calculate remineralization (Lima et al., 2014).To simulate aerosol size distribution dynamics Gelbard et al. (1980) developed a method to solve the Smoluchowski equation they call the sectional method.Jackson and Lochmann (1992) applied the sectional method to marine particles and modeled an algal bloom using 22 size classes to simulate the evolution of the particle-size spectrum of marine aggregates as they coagulate by Brownian motion, shear and differential settling.Burd et al. (2007) expanded this method to study how biological factors and changes in the particle-size spectrum can modify the POC / Th ratio.The sectional method is also utilized in a model of the particle flux in the twilight zone also including settling, microbial activity, zooplankton feeding and fragmentation (Stemmann et al., 2004).Particle-oriented models have found that POC flux is sensitive to aggregation and zooplankton activity (Kriest and Evans, 2000;Stemmann et al., 2004;Gehlen et al., 2006).
We have developed a model, SLAMS, that simulates the flux of orgC and minerals through the water column from the sea surface where production of particles takes place to the seafloor at 4 km depth.We consider two types of particle formation: aggregation of particles (also called marine snow) and grazing and packaging by zooplankton producing fecal pellets.To tune ad hoc parameters, we start with average ocean conditions: SST = 17 • C, PP = 700 mg C m −2 day −1 , bloom index (BI) = 1 and mixed layer depth (MLD) = 50 m, and compare model orgC flux to sediment trap data from Lutz et al. (2002).To see how well the model predicts biogenic flux, we put in respective surface conditions for 11 regions and compare model results to observed orgC flux in the deep, pe-ratio (the fraction of produced organic carbon that is exported) and the rain ratio (orgC / CaCO 3 ) at the seafloor.We are interested in understanding the effects of changing a climatic parameter on the attenuation of the orgC flux.Our model differs from past models in that it has a large stochastic component and the particles are Lagrangian, simulated using a modified super-particle method described in the next section.The advantage gained by the Lagrangian approach is that we can track a large number of aggregate compositions (orgC and minerals) while also dynamically resolving the particle-size spectrum.The goal of this study is to formulate a numerical model intended to simulate these physiobiochemical dynamics.We are not looking to come up with a model to couple to large-scale ocean circulation models but to see if we can piece known mechanisms together and predict fluxes observed in the ocean.Applications for such a model are to tweak parameters and mechanisms to see how the flux responds to environmental conditions to better understand the biological carbon pump.It might also be suitable to develop a parameterization of how the flux responds to external changes in the ecology and environment.Finally, we look at the sensitivity of the model to sea surface temperatures (SSTs), primary production (PP) and a BI, which is a measure of seasonality.2 Model description

Aggregate classes
The constitutive elements of SLAMS are the aggregates.They are clusters of primary particles, a combination of phytoplankton cells (coccolithophorids, diatoms, picoplankton), TEP particles and terrigenous dust particles.SLAMS keeps track of a large number of aggregates that we call aggregate classes (ACs).Each AC carries a suite of information about the state of one class of aggregate in the water column (Table 1).ACs, also called super-particles or super-droplets in atmospheric applications, are representative of some variable number, n, of identical aggregates (Sect.2.5.1).Each aggregate is composed of p primary particles (Fig. 1), forming an aggregate consisting of up to 10 types of orgC (representing different ages), up to 10 types of TEP (also representing different ages) and four types of minerals (calcite, aragonite, biogenic silica and terrigenous material).Lability of orgC and TEP is determined from its age.From this information, SLAMS constructs the physical characteristics of the aggregate that determine its sinking velocity.

Production
In all, 20 new ACs are created per 8 h time step and added to the list of existing ACs.This number is a trade-off between computation time and smoothness of the runs -20 per time step turned out to be more than enough ACs to pro-Figure 1.A schematic of an AC.The complete list of attributes is in Table 1.The AC in this example represents 1000 aggregates each composed of 28 primary particles.
vide consistent results between runs and an acceptable run time.The type of phytoplankton (or terrigenous material) of each new particle is chosen by a Monte Carlo method, wherein a uniformly distributed random number x ∼ [0, 1], is generated and compared with ranges in Table 2 (in this paper all random numbers are uniformly distributed on [0,1]).
Table 2 shows the default functional group in SLAMS that is a typical open-ocean ensemble and would require modification with information on specific functional types and rates of particle production for detailed regional applications.When an AC is produced, p = 1, meaning that the AC represents many copies of a single cell (primary particle) that has not aggregated.For example, if primary production is 1000 mg C m −2 d −1 , each AC (unless it is terrigenous material) initially contains 1000/20 • 8/24 = 16.7 mg C  2) are used for all regions, except the Southern Ocean where calcifiers produce half as much calcium carbonate and diatoms produce twice as much biogenic silica as in the rest of the simulations.
The depth, z, of a new particle is determined by where x is a random number and z is in [cm].Half the production thus takes place in the top 8.2 m.

Bloom index
Blooms are a condition of elevated phytoplankton concentration, possibly a result of complex predator-pray imbalance, and have confounded scientists for decades (Behrenfeld and Boss, 2014).Blooming diatoms have been found to have a lower transfer efficiency than a more carbonate dominated non-blooming flux, suggesting that compositional differences of aggregates and aggregate structure is important for understanding the controls on the organic carbon flux (Lam et al., 2011).Uncertainty in how spring blooms respond to changes in pCO 2 and temperature are still unresolved (Feng et al., 2009).In SLAMS, the BI describes the degree to which production is characterized by blooms: where BI = 1 represents constant production throughout the year and BI = 0.25 means that the annual production all takes place in a quarter of a year.It does not take into effect differences in ecology or physiology of the plankton that may be dominant in blooms vs. not in blooms.Continuing with the example above where PP = 1000 mg C −2 d −1 , if BI = 0.5, then PP = 2000 mg C m −2 d − 1 for the first half of the year and 0 for the second half of the year.

Stickiness and TEP
The kinetics of particle aggregation depend on the encounter rate and the probability of sticking: interparticle attachment rate interparticle collision rate (4) (Alldredge and McGillivary, 1991).The main glue that holds marine snow together in the water column appears to be TEP, a mucus-like polysaccharide material exuded by phytoplankton and bacteria, especially under conditions of nutrient limitation during the senescent phase of phytoplankton blooms (Passow and Alldredge, 1994;Engel, 2000;Logan et al., 1995;Dam and Drapeau, 1995;Kahl et al., 2008).It is exuded in dissolved form but it separates into suspended droplets, or gel, in conditions of turbulence (Verdugo et al., 2004).In SLAMS, TEP production consists of 6 % of the primary production in the default case, in terms of carbon.This number comes from sensitivity studies (Jokulsdottir, 2011) where we simulate equatorial Pacific conditions (SST, PP and seasonality) and compare orgC flux to the deep ocean to sediment trap data.In the model formulation, TEP is not part of primary production but an independent variable that can be changed without altering primary production, allowing for sensitivity analysis of TEP.In each time step, 20 new ACs are produced that are phytoplankton or terrigenous material, and 1 new AC that is TEP.In the example above, 6 % of primary production (1000 mg C m −2 d −1 ) is 60 mg C m −2 d −1 .One TEP particle contains 1 pmol carbon and so n = 5 × 10 9 .The role of TEP in controlling the flux of organic matter is complicated however by its low density (∼ 0.8 g cm −3 ), which acts to decrease the overall density of an aggregate, potentially to the point where it becomes buoyant and ascends rather than sinks (Mari, 2008).In SLAMS, particles rich in TEP can be less dense than sea water resulting in upward movement.Occasional accumulation of TEP-rich aggregates at the sea surface prompted us to consider what happens to buoyant particles at the sea surface and include wind-driven surface mixing in the model.As described above, much of the stickiness of aggregates appears to be due to TEP (Alldredge et al., 1995;Engel, 2000;Verdugo et al., 2004).Our formulation for particle stickiness is simple: TEP particles have α = 1 and at the time of production the stickiness of all other particles is 0. After a particle aggregates, its stickiness is the volume-weighted average of the stickiness of the component particles: where V TEP is the volume of TEP in the aggregate and V a is the volume of the entire aggregate.As a result, in SLAMS, particles do not stick unless TEP is present.

Rate of aggregation
The rate of particle aggregation is governed by two limiting mechanisms: the rate of collision (β) (see discussion in, e.g., Jackson, 1995;Burd andJackson, 1997, 2009) and the rate of sticking (α) once collided.Particles collide by three mechanisms; very small particles mostly encounter each other by Brownian motion, whereas large particles meet most of their partners due to fluid shear and differential settling (i.e., the larger particles settle faster sweeping up the smaller ones).
In SLAMS, the coagulation kernel β is a sum of three mechanisms: collision frequency due to Brownian motion: shear: where and differential settling: (Burd andJackson, 1997, 2009).Here k is the Stefan-Boltzman constant, µ is dynamic viscosity, ν is kinematic viscosity and ε is turbulent dissipation rate set to 1×10 −4 .r i and r j are the radii of the two aggregates being evaluated for collision and u i and u j are the settling velocities of the two aggregates.

Model formulation of aggregation
In the real world, in a given amount of time, some fraction of the particles represented by an AC in our scheme might aggregate with a particular other class of particles and others not.In order to prevent an unmanageable proliferation of particle types, we require that either all of the less-numerous particles find partners in a given time step, or none of them do.The decision is made stochastically using a Monte Carlo method.The idea is that over many possible aggregations, the overall behavior will be statistically similar to a real-world case, or a sectional model, where only a fraction of the particles would aggregate in any given time step.If instead we would allow for a fraction of the aggregates in an AC to aggregate and not all of them, then a new AC would have to be made for every successful aggregation adding hundreds of new ACs each time step.The main advantage of a stochastic over a deterministic approach to aggregation simulation is that here we are able to simulate a large number of compositions without altering the run time.
The number (ξ ) of ACs in a depth bin is variable, depending on aggregation and sinking rate.After the first time step, ξ = 21 in the top depth bin and after the second time step ξ = 42 unless one or more AC sank out.After the model has reached steady state, ξ is a few hundred or a few thousand, depending on the depth and conditions.Each time step, within each depth bin, we pick ξ(ξ − 1)/2 pairs of ACs uniformly at random for potential aggregation.For a given pair, we denote the indices of the ACs by (i, j ), such that n j ≥ n i , i.e., the number of aggregates in AC i is equal or smaller than the number of aggregates in AC j .For each such pair of aggregates, we compute the probability of aggregation between a single i aggregate and n j j type aggregates as P i,j : where α(i, j ) and β(i, j ) are the stickiness parameter and coagulation kernel, respectively.Equation ( 10) represents a continuous time model.In SLAMS we approximate Eq. ( 10) by discretizing in time.As a result, for a given encounter where one of the aggregates is extremely large, it is possible that P i,j > 1.We take a P i,j > 1 as an indication that the time step dt has been chosen to be too large for the approximation to be reasonable hence reducing dt is appropriate.In that case, for the given encounter, the time step is decreased by factors of 10 until P i,j < 1.We derive the expected number of aggregation events between all of the i and j aggregates given by E j : assuming there are n j j particles and n i = 1 i particle and that the probability of aggregation, P , is constant, the expected number of aggregations in time step dt is where q 1 is the number of aggregations the one i particle experiences.For n i = 2, we find the expected number of aggregations that the second i particle undergoes, conditional on q 1 , E q 2 |q 1 = (n j − q 1 )P , www.geosci-model-dev.net/9/1455/2016/Geosci.Model Dev., 9, 1455-1476, 2016 where q 2 is the number of aggregations the second i particle experiences.The total number of aggregations that take place between i and j particles for n i = 2 is For and so the total number of aggregations when n i = 3 is By induction, we obtain the following general formula for the expected number of aggregations when The mean number of j particles that aggregate with a given i particle is therefore Assuming independence of particles, the variance is We round E j to the nearest integer.If E j > 0, aggregates will aggregate.E j j aggregates stick to each i aggregate.
The AC that at the beginning of the time step represented aggregate i now represents n i aggregates with p = p i +E j •p j .The AC that was aggregate j is unchanged after the time step except that it now represents n = n j −E j •n i aggregates.The scheme was designed so that the number of ACs remains unchanged through a aggregation event.This prevents a runaway escalation of the computational load of the run.Figure 2 shows a schematic of two particles sticking.

Physical characteristics of aggregates
An aggregate that forms by aggregation of smaller particles forms a fractal structure with a dimension, D N , that describes its porosity: p ∼ r D N a where p is the number of primary particles and r a is the radius of the aggregate (Logan and Wilkinson, 1990).The closer D N is to 3, the more the structure fills up three-dimensional space, and the lower its porosity.The fractal dimension of marine snow has been inferred from measurements of aggregate properties such as settling velocity, porosity and size, to be in the range of 1.3 to 2.3 (Logan and Wilkinson, 1990;Li and Logan, 1995).The porosity of marine snow appears to be large, always above 0.9 and mostly closer to 0.99, increasing with the diameter of the aggregate (Alldredge and Silver, 1988;Ploug et al., 2008).Fecal pellets are more compact, and thus their porosity is smaller, about 0.43-0.65 (Ploug et al., 2008).The aggregate sinking velocity depends on the aggregate radius and porosity.We follow Logan and Wilkinson (1990) calculating radius and porosity from the fractal dimension and composition of the aggregate (Table 1).The fractal dimension of aggregates in SLAMS is D N = 2.0 (Li and Logan, 1995).Fecal pellets are not fractal in nature, so we use the relationship between fecal pellet volume and sinking rate developed by Small et al. (1979) and find that a porosity of 0.5 results in sinking velocity of model pellets that compare well to that relationship.The resulting model pellet density for 0.5 porosity is near the range of measured pellet density of 1.08-1.2,depending on the composition of the pellet (Ploug et al., 2008;Feinberg and Dam, 1998).The aggregate radius, r a , is obtained from the fractal dimension and the radius of the primary particle, r p : The radius of the primary particle, r p , is calculated from its volume under the assumption that it is a perfect sphere: where V p , the volume of the primary particles in the aggregate, is calculated by dividing the total volume of material in the aggregate, V m , by the number of primary particles: and where X i is the concentration [mol] of substance i, mw i is the molar weight, ρ i is the density and S = orgC, TEP, CaCO 3 , bSi, dust .The porosity, φ, is related to the volume of individual particles that makes up the aggregate to the total volume of the aggregate: where (Logan and Wilkinson, 1990).The density of the aggregate, ρ a , is calculated from the density of the material in the aggregate, ρ s , the density of seawater, ρ sw and the porosity, φ: Note that the total number of primary particles is conserved; n i p i + n j p j = 240 before and after coagulation.Same is true for the amount of organic carbon, CaCO 3 and bSi.

Settling
Sinking velocities of particles range from a few up to hundreds of meters per day (Trull et al., 2008) and have been found to increase with depth (Berelson, 2002).The sinking orgC flux of particles in the ocean, binned according to sinking velocity of the particles, has been found to exhibit a bimodal distribution, with substantial fluxes from particles sinking at about 1 m d −1 , and close to 1000 m d −1 , but very little in between suggesting there are two types of sinking particles that make up the orgC sinking flux (Alonso-Gonzalez et al., 2010;McDonnell and Buesseler, 2010).In SLAMS, the aggregate sinking velocity: is calculated using Stokes' law (Alldredge and Gotschalk, 1988), where g is the gravity of Earth.For low Reynolds numbers where viscous forces are dominant, and where ν is kinematic viscosity.For large Reynolds numbers where turbulence starts to play a role, the drag coefficient is calculated using Whites' approximation: , 1991), which is valid for Re < 2 × 10 4 .Substituting f (Re) with White's approximation results in a nonlinear equation with u as the variable, which we solve for using Newtons' method.Model aggregates only sink vertically, there are no lateral currents.Aggregates reach high Re conditions when settling at a few hundred meters per day, depending on their size.Within the range of salinity found in the ocean, the viscosity of seawater is mainly a function of temperature (Dorsey, 1940).We use a formula that fits empirical data of seawater viscosity at 35 ppt: to find the viscosity, ν, at temperature T [K].The change in viscosity with temperature has a large effect on the sinking velocity of particles.For the change in sea water temperatures from the sea surface to the sea floor in the tropics, the change in viscosity leads to a 50 % decrease in sinking velocity.

Organic carbon and TEP degradation
Bacteria are relevant to the biological pump for their role in degradation of orgC.Organic carbon is a chemically heterogeneous combination of proteins, lignin and cellulose, which vary in structure and degradability (Westrich and Berner, 1984), and as it degrades, its heterogeneity increases.As the structure is altered, it becomes increasingly unrecognizable to enzymes, producing the observed decrease in reactivity T. Jokulsdottir and D. Archer: Sinking aggregates (Emerson and Hedges, 1988).TEP degradation is treated the same way as orgC degradation.Respiration rates have been measured, both on natural aggregates collected in the ocean, and on aggregates formed in laboratory roller tanks from freeze-killed diatom ooze.Higher respiration rates and bacterial production were found in younger (1-3 days) rather than in older (7-14 days) aggregates, as the labile proteins were respired first, leaving the less labile material to be respired more slowly (Grossart and Ploug, 2000).Observed respiration rates of newly produced diatom aggregates range between 0.057 and 0.089 d −1 (Grossart and Ploug, 2001) or 0.083 ± 0.034 d −1 (Ploug and Grossart, 2000), and generally decrease with time, from 0.09 initially to 0.05 d −1 3 days later (Grossart et al., 2003).
To simulate the decrease in degradation rate with orgC age, we track 10 age bins for orgC in each aggregate, and construct a degradation rate that is both a function of age and temperature: where l stands for the age bin and T is temperature The age of the orgC in bin l, age l , is 2 l−2 to 2 l−1 days and This relation assigns freshly produced orgC at the sea surface a degradation lifetime of a few days, and it prescribes a decrease in reaction rate with phytoplankton age (Grossart and Ploug, 2001;Ploug and Grossart, 2000;Grossart et al., 2003) following the observations of Middelburg (1989).The temperature dependence results in about a doubling of the reaction rate with 10 • C of warming, a moderate activation energy of about 45 kJ mol −1 .Aging of orgC and TEP is accomplished by transferring orgC l dt/(age l −age l−1 ) orgC (or TEP) from bin l into bin l + 1 for l = 1, . .., 9. Material does not age out of the top bin, l = 10.If the fraction of TEP in an aggregate is less than 0.02, it breaks into P f fragments as described in Sect.2.9.1.

Zooplankton
The presence of slowly sinking small particles in the deep ocean attests to the effects of disaggregation or fragmentation of aggregates in the water column (Jackson, 1995).In a sectional modeling study, Ruiz (1997) treated disaggregation as driven by turbulent flow in the ocean.Laboratory experiments with diatom aggregates (a fragile form of marine snow) find that stresses in excess of those due to turbulent shear in the water column are required to break them, pointing to biological shear and grazing as the dominant breaking mechanism (Alldredge et al., 1990).The fluid stress required to break aggregates can be found near appendages of swimming zooplankton (Goldthwait et al., 2004).The abundance of marine snow in surface waters appears to undergo daily cycles (Lampitt et al., 1993;Stemmann et al., 2000;Graham et al., 2000), as does concentration of zooplankton.Dilling and Alldredge (2000) observed that the mean size of aggregates decreased when the zooplankton Euphausia Pacifica were abundant.Recently, the idea that zooplankton may break aggregates to decrease their sinking velocity and increase their surface area to encourage bacteria to break down refractory organic carbon and enhance its nutritional value has been proposed as a mechanism to close the carbon budget in the twilight zone (Giering et al., 2014;Mayor et al., 2014).While questions remain unanswered as to the exact effect zooplankton has on marine particles and flux, observations suggest that interaction between zooplankton and marine snow does take place.Zooplankton in our model have the ability to fragment and ingest aggregates and produce a range of fecal pellet sizes.The zooplankton encounter rate is the function in SLAMS that is least supported by science.It is a damped function of how much food there is available.The factor 10 −4.5 is chosen by requiring the orgC export to be 100 mg C m −2 day −2 and the orgC flux to the seafloor to be about 2.5 mg C m −2 day −2 for the default forcing. where is the amount of orgC (food) in each depth bin.For each AC, at each time step, a random number, r, is generated and compared to R enc .If r < R enc , the zooplankton do encounter the AC and either ingest it (all the aggregates it represents) or fragment it.

Fragmentation
To determine whether the aggregates encountered by zooplankton are fragmented or ingested SLAMS first checks if it will be fragmented.We assume the efficiency for breaking small particles is lower than for breaking large ones (Alldredge et al., 1990).We account for this relationship with a parameter, R break , that increases with aggregate radius, so that it is more likely that zooplankton are able to break large aggregates than small, if they do encounter an aggregate.To simulate a smooth function with regard to aggregate size, we came up with the following equation: R break = 0.1 arctan(r a /10 4 ).
If r < R break , the particle fragments into a number of daughter particles.From Goldthwait et al. (2004) we construct a power law with support between 2 and 24 integers (fragments) that describes the number of fragments, f , a particle breaks into: Geosci.Model Dev., 9, 1455-1476, 2016 www.geosci-model-dev.net/9/1455/2016/A random number, r, is compared to P f for each successive number of fragments starting from 2. When r < P f , the number of fragments the aggregate breaks into is found.It is unlikely that very small particles will be fragmented because of how R break is constructed; however, if f > p, then we let f = p meaning the particle fragments into its primary particles.Mass is always conserved during coagulation or fragmentation.When a particle is fragmented, p new = p old /f and n new = n old •f , here p new , n new , p old and n old stand for p and n after and before fragmentation, respectively.If p old cannot be divided into an integer for the given number of fragments, it is divided into the nearest integer below and the remainder mass is divided between the particles and added without adding primary particles.For example, if p old = 17, n old = 1 and f = 3, then p new = 5, n new = 3 and the mass of two primary particles is divided between the three new aggregates.
Here we compromise in conserving the number of primary particles but do conserve mass.

Ingestion
If the particle is fragmented, it cannot be ingested in the same time step, but if it is not fragmented, then it is evaluated for ingestion.To prevent zooplankton from ingesting mineral particles with no orgC and dissolving and repackaging those minerals, an ad hoc parameter is introduced, which adjusts the probability of ingestion as a function of orgC concentration.The orgC weight fraction, orgC wf , is calculated and compared to a random number r.If r < orgC wf , the particle is deemed appetitive and ingested.The model is sensitive to the choice of this parameter.If we choose r < 10orgC wf or r < 0.1orgC wf that decreases the orgC flux by half and increases the orgC flux 3-5-fold, respectively.

OrgC assimilation
The fraction of primary production grazed by microzooplankton estimated using dilution techniques (Landry and Hassett, 1982) leads to 60-75 % of primary production is consumed by protists and between 2 and 10 % is consumed by macro-grazers (Landry and Calbet, 2004).Other studies find that the fraction of production consumed by protists is perhaps a little lower, between 20 and 70 % (Kiorboe, 2000), 22 and 44 % (Riser et al., 2008), and 40 and 60 % (Richardson et al., 2004).The relative rates of carbon consumption in the deep ocean between bacteria and zooplankton vary regionally for reasons that are not well understood (Steinberg et al., 2008;Robinson et al., 2010).By satisfying that about 60-90 % of orgC that is ingested by zooplankton is assimilated and the remaining 10-40 % is egested in fecal pellets (Tande and Slagstad, 1985;Dilling et al., 1998;Bochdansky et al., 1999), the respiration term is both age and temperature dependant: In the model, fecal pellets are represented similarly to aggregates, with the difference that their porosity is specified to be 50 % as described in Sect.2.6.

Biogenic mineral dissolution
pH of the guts of zooplankton increases during feeding, most likely due to dissolution of calcium carbonate (Pond et al., 1995).Jansen and Wolf-Gladrow (2001) modeled zooplankton gut chemistry and conclude that up to 10 % of ingested calcium carbonate can be dissolved during the passage through the gut.Zooplankton mediated CaCO 3 dissolution may help explain the apparent 60-80 % decrease in CaCO 3 sinking flux in the upper 500-1000 ms of the water column (Milliman et al., 1999), and the water column alkalinity source detected by Feely et al. (2004).In SLAMS, dissolution of CaCO 3 in zooplankton guts is related to the fraction of the orgC that is assimilated: Where CaCO 3 (1) stands for calcite, CaCO 3 (2) represents aragonite and kz CaCO 3 (i) is 2 × 10 −2 s −1 for calcite and 4 × 10 −3 s −1 for aragonite.

Abiotic mineral dissolution
In addition to biological dissolution of CaCO 3 in zooplankton guts, there appears to be significant dissolution of the more soluble CaCO 3 mineral aragonite in the water column (Gangstøet al., 2008).A significant fraction of the CaCO 3 production flux is composed of aragonite (Fabry and Deuser, 1991), produced for example by pteropods.Water column conditions reach undersaturation with respect to aragonite at a shallower water depth than for calcite.Biogenic silica (bSi), the other main biomineral in the ocean, is undersaturated throughout the ocean, and is distinguished by a strong temperature sensitivity of its dissolution rate, dissolving significantly faster in warm waters (Gnanadesikan, 1999).In the deep ocean, carbonate dissolution is a function of the degree of undersaturation: where is the saturation state of the sea water for calcite and aragonite.kt (1) = 5 day − 1, kt (2) = 3.2 day −1 , η = 4.5 for calcite and 4.2 for aragonite.The carbonate ion profile or concentration with depth, [CO = 3 ] in situ , is prescribed and dissolution of calcium carbonate does not feed back on it.
For bSi dissolution, the dissolution rate is calculated as and a size-dependent term: to simulate the protective organic matrix or membrane that coats live diatoms and serves as protection to dissolution of frustules (Lewin, 1961;Kamatani, 1982).We assume that non-aggregated diatom cells are alive and thus have a coating to protect them from dissolving.The sensitivity of the flux to k size is investigated in Sect.3.

Sea surface temperature
The temperature is constant in the mixed layer.Below the mixed layer profiles are a linear interpolation from the SST to 4 • C at 1 km to 2 • C at 4 km.In the model, the temperature is constant over the course of the year.

Carbonate ion profile
The carbonate ion concentration is set to 220 µmol kg −1 at the sea surface and 80 µmol kg −1 below 1 km depth.It is linearly interpolated between the surface and 1 km (Anderson and Archer, 2002).

Mixed layer
The mixed layer is isothermal as described in Sect.2.11.Particles in the mixed layer are assigned a random depth within the mixed layer each time step.

A time step
In a time step (Fig. 3) a small number (default = 21) of new particles are produced near the sea surface and their attributes such as composition, sinking velocity and stickiness are set.The water column is divided into depth bins (dz = 10 m).For pairs of ACs within a particular depth range (bin), the probability of aggregation is calculated and compared to a random number to assess whether the pair should stick.If so, the attributes of both particles involved are updated, the more numerous class losing members to aggregation with the less numerous class, which aggregates entirely.The encounter rate of each AC to zooplankton is also calculated for each depth bin.Next, the model loops through every AC and checks if it encountered zooplankton, respires orgC and TEP, dissolves minerals, disintegrates if there is not enough TEP to hold it together, ages orgC and TEP, checks if the size is small and considered to be dissolved, calculates new sinking velocity and settles or ascends accordingly.When the AC reaches the seafloor, its content is cleared and memory is freed for a new aggregate.

Sensitivity studies
To understand the sensitivity of model to parameters not well constrained by experiments, we run the model where a particular parameter is multiplied by a scaling factor (Fig. 4).
For stickiness, α (Eq.5), we look at seven runs where α is multiplied by 1/8, 1/4, 1/2, 1, 2, 4 and 8 and plot the export (Fig. 4a) and flux to the seafloor (Fig. 4b) after the model reaches steady state.Similarly, we do seven runs for rate of bacterial respiration, k(l, t) (Eq.32) and zooplankton fragmentation, R break (Eq.35).The range of TEP production is from 1/6 to 16/6 of the default as the model does not handle well TEP values out of that range.Zooplankton encounter, R enc (Eq.33), is multiplied by factors between 1/8 and 18/8.Stickiness and TEP.The organic carbon flux is very sensitive to the stickiness of TEP and amount of TEP produced (Fig. 4a and b).As α increases (TEP becomes stickier), organic carbon flux is increased.Very little stickiness results in almost a breakdown of the biological pump with very little flux to the seafloor.As stickiness is increased, flux is increased due to an increase in large and fast sinking aggregates.An increase in TEP production works to increase flux up to a certain point, but because of its low density, after that point, more TEP leads to a decrease in the flux.
Bacterial respiration.The fraction of respiration undertaken by bacteria vs. zooplankton is highly variable regionally and between ecosystems (Steinberg et al., 2008;Giering et al., 2014).Here we multiply Eq. ( 32) by a factor between 1/8 and 8.There is very little change in flux until k l,T is 8 times its default value, where it results in an increase in large, fast sinking, fluffy aggregates to the seafloor and a decrease in export (Fig. 4a and b).Zooplankton encounter.Turning off zooplankton entirely results in a scenario where there is almost no flux to the sea floor.In this experiment, particles are relatively buoyant, rich in TEP and sink slowly.When they enter deep waters undersaturated with regard to calcite their density decreases further and they cease to sink.We therefore dialed down the zooplankton encounter rate slowly (Fig. 4c and d).Very little zooplankton encounter results in high export and flux to the sea floor being dominantly large aggregates.As the encounter rate is increased, more small aggregates reach the seafloor and fewer large aggregates, decreasing the flux.
Zooplankton fragmentation.In the model, if an aggregate is encountered by zooplankton, it will either fragment it or ingest it.So, as the probability of fragmenting goes up, less is ingested and respired by zooplankton.The model is not very sensitive to R break .
Diatom organic coating.Figure 5 shows the sensitivity of organic carbon, biogenic silica, calcium carbonate and terrestrial material flux to changes in k size .If diatoms in the model are allowed to readily dissolve (k size = 1), very little bSi makes it to the seafloor and the organic carbon flux is decreased by 50 %.The biogenic silica flux has very little effect on the flux of other minerals.
Particle-size spectrum.The particles size spectrum produced by the model does not seem to be sensitive to environmental or model parameters (Fig. 6).In our model the slope In each panel, the output from a 100 m deep column is plotted.The diagonal lines correspond to slopes −4 (big dash), −3 (small dash) and −2 (dotted).Sensitivity of the slope of the particle-size spectrum at three depths with changing zooplankton encounter rate (P enc ) and stickiness (α) (e) as a function of the scaling factors and slope as a function of TEP production (f).
of the whole spectrum varies around −3 (Fig. 6).At the small end (r < 30 µm) it is less steep (between −2.5 and −3), and at the large end (r > 30 µm) it is steeper (between −3 and −4).This may be attributed to the different mechanisms primarily responsible for coagulation for different size classes (Brownian motion, shear, differential settling) or zooplankton grazing.The slope does not seem to vary systematically to changes in PP, SST or model parameters.
Mixed layer.We investigate the effect of the mixed layer depth on the orgC flux by experimenting with the temperature profile.If the mixed layer is not isothermal and its only role is mixing particles, the flux increases as the depth of the mixed layer is increased.If we let the mixed layer be isothermal and let the thermocline extend from the bottom of the mixed layer to 1000 m, there is no appreciable effect on the orgC flux.Finally, if we let the mixed layer be isothermal and the thickness of the thermocline be constant, 100 or 300 m, the orgC flux decreases as the mixed layer depth increases (Fig. 7).The increased time particles spend in warm waters when the mixed layer is deep contributes to the increased respiration of orgC and thus less flux.

Model-model comparison
The rational for using a Monte Carlo method over a deterministic method is the large number of chemical components that can be resolved without increasing the runtime of the coagulation function.The runtime of the sectional coagulation model with N chemical components can be as high as O(a • b N ), where a and b are the number of size bins and composition-ratio bins; however, with simplifications it is possible to decrease the computational cost (Jackson, 1998).In contrast, for a stochastic model, the runtime of an aggregation simulation is proportional to the number of ACs squared, and does not increase with N.This can be considerable and probably does not make SLAMS suitable for large ocean circulation models without simplifications.
Variations of the Monte Carlo method for aggregation developed by Gillespie (1975) are used to study aggregation in astrophysics (Ormel and Spaans, 2008;Zsom and Dullemond, 2008;Shima et al., 2009), atmospheric chemistry (Riemer et al., 2009;Shima et al., 2009) and oceanography (Khelifa and Hill, 2006;El Saadi and Bah, 2007).The approaches vary to fit the goal of each study, for example, in  (Wetherill, 1990) (solid lines) at times 1, 2, 4, . . ., 64 and model results for corresponding times: 1 (plus-signs), 2 (crosses), 4 (asterisk), 8 (open boxes), 16 (solid boxes), 32 (open circles), 64 (solid circles) (a).In all, 10 000 ACs were used here and a total of 1 × 10 6 particles.Results were averaged over 50 runs to generate this plot.For this comparison, particles do not sink and the collision kernel is a constant (1 × 10 −6 ).The vertical axis is the number of particles of mass M at time t times M 2 .Variation of total mass with time (b) with the collision kernel used in SLAMS for eight different initial concentrations.Particles coagulate and sink but are not altered by biological or chemical processes.
whether mass or number of ACs is conserved, and in how many true particles an AC represents.
One computational difference between the our model and some of these is that we use a super-droplet method (Shima et al., 2009;Ormel and Spaans, 2008;Zsom and Dullemond, 2008), where each simulated particle (AC) represents a large and variable number of real particles.The super-droplet method bypasses computational problems that have historically been considered a great drawback of Monte Carlo methods.Another distinction is that our model allows for multi-stage aggregation within a given time step.The main focus of our model of sinking aggregates is tracking particles, their composition and geometry, as they aggregate and are chemically and biologically altered by bacteria, zooplankton and water chemistry.
To assess the validity of the aggregation method derived in Sect.2.5.2, we first look at the evolution of the particle spectrum with time where particles are not allowed to sink and the collision kernel is set to constant.We compare this spectrum to the analytical solution as presented by Wetherill (1990)  (Fig. 8a).The Monte Carlo algorithm is able to reproduce the time evolution of the particle-size spectrum when a minimum of 100 aggregate classes are tracked.(For the full model runs described below, the number of computational particles reached about 30 000.)As the statistical significance of the larger particle-size class in the Monte Carlo model declines, its abundance is subject to random fluctuations, as seen in the large variations in the Monte Carlo model at t = 64 (Fig. 8a).
In the second test we allow the particles to sink and the curvilinear collision kernel used in SLAMS is used to determine aggregation.We look at the decrease of the mass remaining in the top box with time as also studied by Burd and Jackson (1997) (Fig. 8b).The similarity of the simulation results boosts our confidence in the Lagrangian model.

Model-data comparison
Particle-size distribution is a property of marine system; it provides information about the structure of the ecosystem and particle dynamics.The particle-size spectrum in the ocean has been measured using particle counters and imaging methods.The slope of the particle-size spectrum is usually calculated by dividing the concentration C of particles in a given size range r by the size range: n = C/ r.A global synthesis of the particle-size spectrum found its slope to be in the range −2 to −5 (Guidi et al., 2009) and generally that low PP regions were associated with a steep spectrum, whereas greater productivity regions often exhibits a slope between −2 and −4.In a study confined to the Arabian Sea, the slope varied between −2 and −4 and again the greater negative values were generally found where PP was lower (Roullier et al., 2014).In our model the slope of the whole spectrum varies around −3 (Fig. 6).
To see how well the model captures the flux of material in the water column, we compare three model results to data: the orgC flux in the water column, the pe-ratio and the rain ratio at the sea floor.Sediment trap data exist in low resolution, temporally and spatially.Existing data have been analyzed quite extensively (Lutz et al., 2002;Dunne et al., 2005;Honjo et al., 2008).Our comparison of data with the model is based on the synthesis of Lutz et al. (2002), who looked for correlations between primary production and export flux and the attenuation coefficient b in the Martin curve (Eq. 1) of trap fluxes from 11 different regions of the world's oceans.Table 3 shows the values of model driving parameters for the 11 regions and Fig. 9 compares model results with data.Ex-   2. cept for the Panama Basin, the model predicts the flux to the seafloor within a factor of 2-5 (Fig. 9).The topography of the Panama Basin is complex and lateral transport of resuspended material has been observed, possibly explaining the high orgC flux seen in sediment traps (Bishop et al., 1986).
In SLAMS, export is defined as flux out of the top 100 m.To see if the model captures relationships to do with export seen in the real world, model responses to sea surface temperature and rate of primary production are compared with data in Fig. 10.These data were compiled from field observations of primary production and particle export in Dunne et al. (2005).SLAMS produces the relationship between SST and export production quite well (Fig. 10a).The positive relationship seen in the date between export production and high PP is also produced by the model, but because in most of our test cases PP is relatively low this relationship is not demonstrated very clearly in Fig. 10d.No relationship is seen in the data between export production and SST, as well as between pe-ratio and PP.The model also sees no relationships between these pairs (Fig. 10b and c

respectively).
Lastly, we compare the rain ratio (organic : inorgC) in the sinking carbon flux at the seafloor with rain ratio data (Honjo  et al., 2008).SLAMS reproduces the observed value of approximately 1 in the tropics with higher values in the subtropics and polar latitudes (Fig. 11).The organic : inorgC ratio of primary production is the same in all model runs except for the Southern Ocean, run where biS production is doubled but CaCO 3 production is halved.the fraction that reaches the seafloor (the seafloor ratio) as a function of these parameters.The model responds to all parameters.An increase of 1 • C results in 1-1.6 % decrease in the pe-ratio (Fig. 12a).
The sea floor is also sensitive to changes in SST: as SSTs drop more orgC makes it to the sea floor.It is most sensitive at low SSTs.Where SST is lower than 10 • C large (> 1 cm) and fast (> 100 m d −1 ) sinking particles contribute significantly to the flux (Fig. 13a and b).Given other things are equal, a lower SST increases the efficiency of the biological pump.
An increase in primary production increases the number density of particles and the rate of aggregation.More large particles are produced that have high sinking velocity where PP is high (Fig. 13c and d).OrgC is exported and transported to the seafloor more efficiently where PP is high (Fig. 12c  and d).For example, comparing PP at 1400 with 1800 (29 % increase in PP) results in a 97 % increase in orgC flux to the sea floor, primarily contributed by particles at the large end of the spectrum (Fig. 13c and d).
Both the pe-ratio and the seafloor ratio are sensitive to the bloom index.Blooms transport orgC more efficiently to depths.Organic carbon flux to the seafloor is approximately quadrupled when annual production takes place in one-third of a year compared to when it is spread evenly throughout the year.The effect of concentrating production in blooms is the same as increasing PP; an increase in cell number densities enhances aggregation and produces more large and rapidly sinking aggregates (Fig. 13e and f).
Our results suggest that the biological pump is most efficient in environments that have low SSTs and high or sporadic primary production, such as polar environments.Our model shows the biological pump is also most sensitive in polar conditions.For example, in a global warming scenario where an increase in stratification leads to a smaller nutrient supply to the surface and a decrease in PP, it predicts greater organic carbon transport decrease in a polar environment than in a tropical environment.Similarly, a 2 • increase in SST would result in a greater decrease of orgC transport in a polar environment than in a tropical environment.

Conclusions
We present a new computational approach to simulating the biological pump from the surface ocean to the deep sea by resolving collections of individual particles.The method allows the model to resolve detailed attributes of the particles that may affect their sinking/degradation dynamics, includ-T.Jokulsdottir and D. Archer: Sinking aggregates ing the mineral fraction, the age distribution of the organic matter and the impact of transparent exopolymer particles (TEP) on particle aggregation.The model is able to reproduce sinking fluxes of orgC through the water column observed in sediment traps, and some of its regional variation.The model fluxes are sensitive to the three driving parameters analyzed: SST, primary production and bloom index.The mechanistic links between these climate-sensitive drivers and the physics and chemistry of sinking particles in the ocean will allow us a better understanding of the direction and magnitude of an ocean biological pump carbon cycle feedback to climate change, both in the past and in the future.
We are exploring the feasibility of using our methodology in a larger model.A study is under way to simplify SLAMS significantly, but keep the Lagrangian aspect, and couple it to a three-dimensional ocean circulation model.It would be interesting to look at the importance of the depth of the euphotic zone on the flux and compare to the findings of Buesseler and Boyd (2009).Transport efficiency is a metric that has been used to describe the variability of the biological pump (Lam et al., 2011;Henson et al., 2012;Marsay et al., 2015) that would be interesting to explore with SLAMS.A future step in the development of SLAMS is to add nutrientand light-limited production.We would also like to put in a more mechanistic representation of zooplankton activity and see if we can better understand the effect of the relative importance of microbial vs. zooplankton respiration.We should like to do more detailed regional studies and perhaps parameterize b in terms of environmental variables such as SST, PP, depth of the mixed layer and euphotic zone or ecology.

Figure 2 .
Figure 2.An example of two aggregate classes sticking.Top left shows ACs i and j before coagulation takes place.Top right a schematic of the aggregates the AC represent.Bottom three rows shows ACs if E j = 1, 2 or 3. Note that the total number of primary particles is conserved; n i p i + n j p j = 240 before and after coagulation.Same is true for the amount of organic carbon, CaCO 3 and bSi.

Figure 3 .
Figure 3.An outline of what the model does in a time step.In parenthesis are the sections where each task is explained.

Figure 6 .
Figure 6.Model output of particle number spectrum at three depths: 100 m (a), 1 km (b) and 4 km (c).In each panel, the output from a 100 m deep column is plotted.The diagonal lines correspond to slopes −4 (big dash), −3 (small dash) and −2 (dotted).Sensitivity of the slope of the particle-size spectrum at three depths with changing zooplankton encounter rate (P enc ) and stickiness (α) (e) as a function of the scaling factors and slope as a function of TEP production (f).

Figure 7 .Figure 8 .
Figure 7. Temperature profile for the seven mixed layer depths investigated (a).Organic carbon export (b) and flux to the seafloor (c) as a function of mixed layer depth (MLD).Organic carbon flux down the water column for the seven depths investigated (d).

Figure 9 .
Figure 9.The flux of organic carbon in 11 regions.This model (solid line) and data (symbols) from Lutz et al. (2002) and two Martin curves using b estimated from Guidi et al. (2015) (see Table3) (line-dot) and b = −0.858(dotted).The model reaches steady state at about 5-10 years.We run two instances of each experiment for 18 years and take the mean of the last 4 years of the two runs.

Figure 10 .
Figure 10.Pe-ratio (the fraction of produced organic carbon that is exported) as a function of SSTs (a) and as a function of PP (b) and export production (orgC flux out of 100 m) as a function of SSTs (c) and PP (d) in tropics (plus signs), subtropics (crosses), subpolar (stars) and coastal (open squares) regions.Data fromDunne et al. (2005).Model results (solid squares) from the 11 regions listed in Table2.

Figure 11 .Figure 12 .
Figure11.The rain ratio (organic carbon to inorganic carbon) of material reaching the seafloor as a function of latitude.Data from the Southern Ocean (circles), Pacific Ocean (open squares), Indian Ocean (stars), Atlantic Ocean (crosses) and Arctic Ocean (plus sign) (data fromHonjo et al., 2008).Model results (solid squares) for the 11 regions mentioned earlier.

Table 1 .
Attributes of a particle that are kept track of by the aggregate class.

Table 2 .
The probability, P , for a type of particle to be produced and range of values compared to a random number to determine the type produced, the amount of material in each particle [pmol particle −1 ], n is the number of particles produced each time.Radius, r [µm], density, ρ [g cm−3].(Particle types: C is coccolithophorid, A is Aragonite forming phytoplankton, D is diatom, Pi is picoplankton, D = dust.