A discrete interaction numerical model for coagulation and fragmentation of marine detritic particulate matter (Coagfrag v.1)

. A simpliﬁed model, representing the dynamics of marine organic particles in a given size range experiencing coagulation and fragmentation reactions, is developed. The framework is based on a discrete size spectrum on which reactions act to exchange properties between different particle sizes. The reactions are prescribed according to triplet interactions. Coagulation combines two particle sizes to yield a third one, while fragmentation breaks a given particle size into two (i


15
The biological carbon pump is responsible for a significant fraction of the organic carbon exports from the surface to the deep ocean (Passow and Carlson, 2012;Le Moigne, 2019), thereby influencing the climate (Kiørboe and Thygesen, 2001). Questions regarding the quantification and prediction of its efficiency and response times are still broadly unanswered. The carbon pump obeys to a wide variety of processes, involving the co-action of a large number of physical, chemical and biological variables (Denman et al., 2007). Coupled Ocean General Circulation Models (OGCMs) and Biogeochemical models (BGCs) contribute 20 to our understanding of the relative importance of these processes. They were developed for decades to assess the biological pump dynamics at the global scale (e.g. Palmer and Totterdell, 2001;Aumont et al., 2003;Dutkiewicz et al., 2005) or for specific regions of the world (e.g. Sarmiento et al., 1993;Blackford and Radford, 1995;Doney et al., 2002;Wiggert et al., 2006;Karakaş et al., 2009).
Since the pioneering work of Riley (1946), BGCs have been widely used in oceanography and their complexity never 25 ceased to increase (Vichi et al., 2007;Anderson and Gentleman, 2012). As reported by Leles et al. (2016), they evolved from Nutrient-Phytoplankton-Zooplankton-Detritus type (NPZD) (e.g. Palmer and Totterdell, 2001), where multiple species and types of organic and inorganic constituents of the pelagic system are gathered into a broadly defined trophic compartment (Doney et al., 2003), to food-webs of multiple Plankton Functional Types (PFTs) (Hood et al., 2006;Anderson et al., 2010).
In most coupled OGCMs-BGCs, efforts have concentrated on achieving an accurate representation of primary and secondary 30 producers-related particle dynamics, while the detritic compartment is generally reduced to just one variable. Detritus are nevertheless essential in the downward export of carbon (Hill, 1992;Kriest, 2002). Considering only one detritus variable may lead to important biases in carbon flux estimations. The size diversity of marine particles is wide, ranging from large, rapidly sinking particulate material (i.e. marine snow) to small suspended particles and relatively non-labile dissolved organic matter (i.e. colloids), which usually also show the highest abundances. Then, their representation through a unique variable and 35 mean values of its descriptive parameters such as sinking velocity (Doney et al., 1996;Lima et al., 2002;Aumont et al., 2003;Dutkiewicz et al., 2005;Kishi et al., 2007) to cover the entire size range in BGCs is questionable.
To overcome this caveat and improve carbon export assessments, the number of detritus-related variables in BGCs models can be increased to better represent the diversity of the sinking particulate matter. For example, some studies such as Moore et al. (2002); Wiggert et al. (2006); Yool et al. (2011); Butenschön et al. (2016); Kearney et al. (2020) defined multiple (two 40 or more) detritic compartments, each being connected differently with other variables and having constant settling velocities.
This approach can certainly increase the level of realism, as these parameterizations are made based on field or experimental evidences (Doney et al., 2003), but we see two major problems with them. First, the description in the numerical framework of a high number of state variables enhances BGCs models complexity, and augments the number of parameters required to characterize relations among those variables (Denman, 2003). Ultimately, it makes it challenging to properly couple them to in aggregate formation. Fragmentation is the opposite process, breaking particles into smaller pieces. Both processes are then affecting particle size distribution but are barely explicitly implemented or parameterized in OGCMs-BGCs.
It is indeed based on the seminal work of Gelbard et al. (1980) on the sectional representation of aerosol size distribution evolution due to collision and coagulation events, that the first coagulation models applied to marine snow emerged. Jackson 60 and Burd (1998) extended the model of Gelbard et al. (1980) and applied it to the marine environment, pointing out the role of fragmentation to counter balance the importance of coagulation (Jackson et al., 1995;Hill, 1996). These innovative works have been widely used in various ecological contexts and linked to plankton models as listed by De La Rocha and Passow (2007) and Jackson and Burd (2015) (e.g. for coagulation refer to Jackson, 1990;Kriest and Evans, 1999;Jackson, 2001;Kriest, 2002, and for fragmentation to Alldredge et al., 1990;Dilling and Alldredge, 2000;Kiørboe, 2000;Ploug and Grossart, 2000;Goldthwait 65 et al., 2004;Stemmann et al., 2004) but scarcely coupled to OGCMs. It is probably due to the remaining quantitative unknowns regarding the ensemble of processes affecting transport efficiency of the particulate organic matter to depth by constraining particle size distribution even after decades of extensive work (Le Quéré et al., 2005;De La Rocha and Passow, 2007). The intrinsic heterogeneous nature of these processes at all spatiotemporal scales increases the challenge to properly implement them in complex OGCMs-BGCs and to evaluate and predict the ocean's role in the Earth's carbon budget. Similar challenges 70 arise in atmospheric GCMs. According to Kang et al. (2019) GCMs with full cloud microphysics are still at an early stage in terms of understanding and simulating many observed aspects of weather and climate, and research is needed to circumvent these difficulties.
In order to circumvent the issues related to the representation of the dynamics of the complete particles size spectrum for OGCMs, we develop in this study a new numerical framework where a particles' size range is discretized in size bins, and where 75 concentration dynamics over these bins is driven by coagulation and fragmentation reactions. This framework is designed to conserve mass over the size range and reactions, and can accommodate size and mass linear and nonlinear discretizations.
Since the main motivation for developing such a model is to provide a tool allowing to characterise detritic variables relations and dynamics in coupled OCGMs without unreasonably increasing the computational cost, a formulation is sought that will attenuate the dependence of the results to the size discretization resolution. Numerical experiments are designed to study the 80 dependence of the results on i) the number of size bins used to discretize a given size range (i.e. the resolution) and ii) the type of discretization (i.e. linear vs nonlinear). Innovations of the approach as regards to previously developed coagulationfragmentation models and designed detritic state variables is briefly discussed as well as the potential for further improvements to allow the inclusion of the presented model into OGCMs to better estimate current carbon export.
The outline of the paper is as follows. Section 2 presents the model description, Sect. 3 describes numerical experiments and 85 Sect. 4 presents and briefly discusses the results. A summary of our main conclusions is given in Sect. 5.

Discrete size spectrum
Whereas most laboratory and field studies estimate particle' number concentration, n ( m −3 ), along a size spectrum (McCave, 1984;Jackson et al., 1995), OGCMs use tracers (or element') concentration, C (mmol m −3 ), (e.g. carbon or nitrogen) to study 90 fluxes among the model compartments (Doney et al., 1996). These two variables are linked by where N is the particle' content of the chosen currency (mmol). Considering a closed system in a given water volume with no particle sources or sinks, particles may be sorted over a size range L s with s a chosen size property (e.g. diameter or volume). To transpose the variables from the continuous form Eq. (1) 95 to a discrete form, L s must be discretized in size bins p such that where ∆s p = sp+1 sp ds is the size range of the bin p, and (s p , s p+1 ) are the lower and upper size bounds of this bin. Therefore, the particle' content of a given bin p can be interpreted as its mean value

100
Note that other particle' properties (e.g. diameter, volume, density, etc.) are similarly interpreted as a bin-averaged value. For example, the particle' diameter corresponding to a given bin can be defined by D p = 1 ∆sp sp+1 sp Dds. This diameter can in turn be used to determine the particle' volume using an allometric relationship such as, V p = λ 1 D λ2 p , (Jackson et al., 1997;Li et al., 1998;Zahnow et al., 2011) 1 . Using this bin-averaged volume (V p ), the particle' content can be redefined as where λ 3 and λ 4 are parameters empirically determined through field and laboratory experiments depending on the element chosen (Alldredge, 1998, see Table 1).
The discrete form of Eq. (1) thus becomes where n p is the particle' number concentration in p (i.e. the number of particles in p), and N the total number of resolved 110 size bin. In a closed system without sources and sinks of particles, the total concentration, C T = N p=1 C p , is required to be conserved over time. The time evolution of the concentration inside a given bin obeys the simple differential form where R p represents all reactions occurring in p.

115
Coagulation and fragmentation are, by essence, reactions that involve three particles. Coagulation involves two particles, with indices i and j, that collide and stick together to form a third, larger one with index k. Conversely, fragmentation can be considered as the opposite reaction where a particle k breaks into two smaller ones i and j, in line with what was observed by Alldredge et al. (1990) in laboratory experiments. Reactions involving more than three particles can always be decomposed as a sequence of triplet reactions. Note that i and j can originate from identical or different size bins. In a linear size discretization, 120 by definition, the k index always refers to a particle belonging to a different bin than particles i and j. Then, the volume V k of particle k resulting from the coagulation of particles i and j = 1...k obeys to: The reaction is arbitrarily built so that the i th particle is always the smallest one. From this assumption, the coagulation reaction can be written as (Fig. 1). However, in the case of a nonlinear size 125 discretization, this rule might be violated and this situation is discussed in Sect. 2.6.
In this model, the rules described above for particles will be applied on the bins discrete spectrum of Eq. (5). For example, in a coagulation reaction implying a given triplet of bins (i, j, k), the rate of change of the concentration in i will depend on that of j, and they will conjointly prescribe the rate of change of k. Such a reaction for a triplet of bins can be interpreted as 130 multiple reactions for triplets of particles, that can be found in their respective bins (i, j, k). For a reaction implying a given triplet of bins (i, j, k), the evolution of the concentration in each bins is described by the following set of differential equations: δC k i,j is a triplet operator that represents both coagulation and fragmentation reactions acting on a given bin :

135
A bold index indicates the bin on which the reaction applies (see also Fig. 2 for a visual explanation of this convention). By construction, of the total number of k particles involved in the reaction in Eq. (8), half of this number is associated with i and the other half with j ( Fig. 2), which explains why we multiply all the terms in Eq. (9) by 1/2. K is the coagulation rate, while F is the fragmentation rate. Coagulation has been studied extensively in previous works both in atmospheric (Pruppacher and Coagulation is the process by which two particles with indices i and j collide and stick together to form a third, larger one with index k. Fragmentation is the opposite process by which a particle k breaks into two smaller ones i and j. Size bounds are shown by vertical dashed gray lines. Reactions involving two small particles (a) from the same size bin (i = j) and (b) from two different size bins (i = j) are shown.  Jackson et al. (1997) [b] Alldredge (1998) [c] Li et al. (2004); McCave (1984) three main collision mechanisms (Kiørboe et al., 1990;Ackleh, 1997;Engel, 2000;Jackson, 2001): Brownian motion, fluid velocity shear and differential settling. Fragmentation (F) can be driven by biology (e.g. related to zooplankton activities such as grazing (Banse, 1990;Green and Dagg, 1997), and swimming behavior (Dilling and Alldredge, 2000;Stemmann et al., 2000;Goldthwait et al., 2004)) or driven by physics (e.g. scales of turbulence (Alldredge et al., 1990;Kobayashi et al., 1999)).

Reaction Matrices
Let us now consider a set of discrete bins (p) that are linearly incremented and have indices (i, j) running from 1 to N . By construction, this yields reactions for k ranging from k = 2 to 2N (i.e. 1 + 1 to N + N ). To account for the concentration evolution associated with all possible reactions, we define four matrices built from the triplet operator (δC k i,j ) defined earlier, D i=j , B i,j , T j,i and F i,k , referring to Diagonal, Bottom, Top and Final matrices, respectively.
Matrix D i=j has dimensions N × N and accounts for reactions in which the two particles have the same size (i = j) The square brackets above indicate the matrix dimensions. Matrix B i,j accounts instead for all reactions acting on i only while T j,i accounts for those acting on j only (where j > i, Fig. 2b,d)

160
In F , the double line separates the resolved size range (above) from the unresolved one (below). The latter contains reactions involving particle sizes outside the resolved size range of the spectrum (i.e. reactions for which k > N ).

The unresolved range
Solution strategies to parameterize reactions in the unresolved range must obey two basic rules: i) they must conserve the total concentration (C T ) in absence of external sources and sinks of particles and ii) they must limit the unbounded growth 165 of the particle size due to coagulation. Here we propose a simple closure to account for the reactions in this range. In order to comply with conservation of the total concentration in the size range (L s ), at least one new bin must be added in which concentration fluxes in and out of the resolved range are stored. Moreover, in order to avoid unbounded growth of the size range due to coagulation, this extra bin must not be allowed to further coagulate with itself or any of the other particles sizes.
As such, all reactions that fall into the unresolved range will be accounted for in a single additional bin for which coagulation 170 is prohibited but fragmentation back into the resolved range is allowed. This extra bin can thus be interpreted as an average of all the reactions in the unresolved range (N < k ≤ 2N ) and will be referred to as k = 3 2 N . Applying this to Eqs. (10), (11),

175
where ∆C = 2δC when i = j and δC otherwise, and K = 0 in ∆C for k = 3 2 N . Equation (14) includes all reactions acting on either i or j, while Eq. (15) includes all reactions acting on k. Since the unresolved range only involves reactions acting on k, all the terms below the double line in Eq. (14) are set to zero. Moreover, each term of the last row in Eq. (15) can be viewed as a sum over all the elements of a given column below the double line in Eq. (13). Note that, for simplicity, we choose to add only one extra bin in the unresolved range, but one could alternatively choose to add up to N bins in order to improve this 180 parametrization.

Summary of all reactions
Based on matrices (14), (15) and specific reaction rules previously described, the reaction vector R p , representing all the reactions for a given bin (1 < p < N + 1) is given by Coag and Frag in Eq. (16) are used to explicitly show the sign of the coagulation and fragmentation reactions from Eq. (9). In other words, coagulation is removing concentration from p in X p + Y p , while adding concentration to p in Z p (and conversely for fragmentation). Equation (16) can be rewritten as a sum of the following series: where ∆C = 2δC when p = 2i and δC otherwise. An example of a complete set of reactions with N = 4 and its additional bin 3 2 N = 6 is shown in Fig. 3. Focusing on the resolved range only, Eqs. (16) and (18) can be combined to yield an alternative formulation of the discrete Smoluchowski equation (Smoluchowski, 1916) 195 where δ ip is the Kronecker delta function that is equal to 1 for i = p and zero otherwise. Notice that the above equation gives the rate of change of concentration, whereas the traditional formulation for the Smoluchowski equation is written in terms of number of particles. Equation 19 can thus be reformulated in terms of number of particles as: 200 The factor 1/2 in the second term ensures that the combination of two particles yields a single larger particle. This is in contrast to the concentration, for which the combination is additive (see Eq. 8).

Robustness to resolution
For a discrete representation of the full size range, L s , larger values of N imply higher resolution of the size spectrum. The total number of reactions of the system described above increases with the size of the reaction matrices, N (N + 1), which is nearly a quadratic function of resolution. On the other hand, coagulation and fragmentation are respectively quadratic and linear functions of concentration, as shown in Eq. (9). Since concentration is itself a quantity that depends on resolution (Eq. 5), an asymmetric response between coagulation and fragmentation to changes in resolution is expected.
In order to build an intuition on the effect of resolution on the reactions, consider a conservative system with a given total concentration, C T = N p=1 C p , and a linear size discretization such that both the particle' content and concentration 210 are constants given by N p = 1 and C p = C T /N , respectively. To further simplify, assume that the rate of coagulation and fragmentation are constants given by K and F respectively. For a conservative system, the sum of all reactions integrates to zero, i.e. N +1 p=1 R p = 0, such that C T remains constant at all time. Since the sum of all reactions does not allow to keep track of the total concentration exchanged along the spectrum, we instead define from Eq. (16) the total reaction amplitude as 215 and since N = 1, where the first term on the right hand side is the total coagulation amplitude and the second term is the total fragmentation amplitude. In this simple example, dependence to resolution is revealed through the respective prefactors N +1 N and (N + 1).
It thus becomes obvious that fragmentation is more sensitive to resolution than coagulation. This result can be explained a 220 posteriori if we notice that the total number of reactions is nearly quadratic with N , while the coagulation concentration is proportional to N −2 , cancelling most of the variation with N . In contrast, the fragmentation concentration is proportional to N −1 , which yields a larger residual dependence on N . To counterbalance this dependence on resolution, the reaction terms are divided by their respective resolution-dependent coefficients such that For simplicity, we chose to describe the above framework using a linear size discretization (i.e. where bins are equally distributed along the size range). However, this choice was arbitrary and we now generalize the framework to a nonlinear size discretization, which is a more natural choice for representing marine particles. A nonlinear size discretization can be seen as local variations of the resolution in the full size range. For example, for a given total number of bins, N , switching from a linear to a logarithmic spacing increases resolution for the small particles, while decreasing it for the large particles. Intuitively, this 235 choice seems better suited to represent marine particles that have a wider variety of microscopic than macroscopic particles.
In this context, the main difference with the framework described in the previous sections is that volume conservation can be violated when using a nonlinear size spectrum, i.e. Eq. (7) is no longer valid. Although counter-intuitive, this does not imply however that mass conservation is necessarily violated. Consider, for example, a simple nonlinear discretization where bins are each separated by an order of magnitude (1 µm 3 , 10 µm 3 , 100 µm 3 , etc.). Coagulation of particles belonging to bins 1 µm 3 and 240 10 µm 3 would ideally produce particle size 11 µm 3 . However, since 11 µm 3 is much closer to bin 10 µm 3 than bin 100 µm 3 (the next larger bin) all the concentration associated with this reaction will fall into the 10 µm 3 bin, thus violating volume conservation, yet conserving the concentration associated with the reaction. We thus modify Eq. (7) to allow that the resulting size of a coagulation reaction is not required to be strictly equal to the sum of the two reacting particles (and conversely for fragmentation), i.e.
The main consequence of Eq. (23) is that it modifies the reaction matrices (14) and (15). In Eq. (14), only the k indices will be modified by a nonlinear discretization. In Eq. (15), the elements themselves will be redistributed on different rows of the matrix. For a logarithmic discretization that enhances resolution towards smaller particles, elements will be moved upwards in the F matrix as an increased number of reactions yield V i + V j = V i+j . For example, δC 3 1,2 could be switched from the 250 third row using a linear spectrum to δC 2 1,2 in the second row using a logarithmic spectrum. Conversely, elements will be moved downwards in the F matrix if a discretization that enhances resolution towards larger particles is chosen. Here, we do not explicitly write the matrix encompassing all possible cases since this would be unnecessarily complex. The construction of the reaction matrices (14) and (15) is done numerically at the beginning of our algorithm (see Gremion and Nadeau (2021)).

Application to a plankton ecosystem model 255
The framework presented here gives rise to N variables that together represent detritic particulate matter. Like other variables of plankton ecosystem models, these additional variables must be treated as Eulerian tracers submitted to diffusive and advective transport. In a typical three-dimensional OGCM, the evolution of the carbon concentration C p (mmolC m −3 ) belonging to the size bin p, is given by where u =îu +ĵv +kw is the velocity field and ∇ =î ∂ ∂x +ĵ ∂ ∂y +k ∂ ∂z is the nabla operator. The terms on the right side are: (i) the advective flux convergence, (ii) the diffusive flux divergence, with K the turbulent diffusivity, (iii) the background vertical advection due to the settling velocity w p associated to size bin p, (iv) the sources and (v) the losses of detritic matter associated with other biogeochemical processes, and (vi) the reaction term, R p , representing coagulation and fragmentation derived in this paper (Eq. 16). The three dimensional velocity u is provided by equations driving geophysical fluid dynamics. The vertical 265 settling velocity, w p , is here assumed constant for a given size, but can vary considerably from one size to another as it strongly depends on particle' properties such as its volume, density and porosity. Therefore, C p does not vary due to differential settling (divergence or convergence), but C T can through the action of the reaction term, R p . In order to focus uniquely on the reaction term, R p , we do not solve Eq. (24) explicitly in this work and leave this for a subsequent study. In the following we use the simplified zero dimensional form ∂C p /∂t = R p to investigate the robustness of the proposed framework on the resolution, N .

270
3 Numerical experiments Two model configurations are set up to study detritic carbon concentration (C) dynamics experiencing coagulation and fragmentation reactions using the previously described model. The first configuration, named E1, uses a linear size discretization over the arbitrarily chosen range 0 to 8 (Fig. 4a). The second configuration, identified as E2, uses a nonlinear size discretization that more realistically represents particle' number distributions observed in the ocean (i.e. particle' size from D 1 = 1 µm to 275 D N = 1 cm, (Stemmann et al., 2004;Monroy et al., 2017), Fig. 4b). Each set of experiments within the configurations, is performed using two different numbers of size bins (N ), in order to study the impact of resolution on the model. The high res-value of 3 2 N is added to represent the unresolved size range, which increases the total number of bins to N = 401 and N = 5, respectively. For each resolution, three simulations are performed: two simulations where coagulation (K) and fragmentation 280 (F) are considered separately, and one simulation where they are combined (KF). In total, twelve simulations are performed with parameter values summarized in Tables 1 and 2. The low resolution (LR), in red, has N = 4 size bins plus one unresolved size bin 3 2 N = 6, while the high resolution (HR), in green, has N = 12 size bins plus one unresolved bin 3 2 N = 18. This is an illustrative example as the HR simulations performed in this paper with N = 400. Size bin bounds are shown by vertical dashed gray lines. All simulations are initialized with the same total carbon concentration CT (see Table. 1). The initial concentration is spread uniformly over the resolved range for E1-LR and E1-HR, and following a power law for E2-LR and E2-HR (see Table. 2). No concentration is initialized in the unresolved 3 2 N size bin. As HR has more size bins than the LR, concentrations values are consequently lower as determined by Cp = C T N in E1. Concentrations are prescribed to the middle size value of a given bin.

Initialisation
For each configurations, all experiments are initialized from a reference distribution at very high resolution (i.e. N † = 4000 bins) that is meant to represent an ideal distribution. This discretization is indicated by the symbol †. The total carbon concen-285 tration is set to C T = 100 mmolC m −3 , and initial concentrations, C p (t = 0), are then initialized for HR and LR by integrating on the reference distribution using a discrete version of Eq. (3) where p refers to indices of the HR and LR discretizations and q to indices of the reference distribution. For HR and LR in configuration E1, the initial carbon concentrations of the reference is uniform (i.e. C † p = C T N † ) as for the particle carbon content 290 (N p = 1). The chosen time step is ∆t = 86400 s and it is run for one time step.
In E2, the reference is initialized with a power law distributed carbon concentration of the form (McCave, 1984;Li et al., 2004) where n is the number of particles of diameter D and Ψ is the slope of the distribution, by connecting with Eq. (5). Unlike 295 in E1, particle' carbon content is set to a size-dependent function following Eq. (4). HR and LR initial carbon concentration distribution are then obtain following Eq. (25), and the model is integrated for a day with a timestep ∆t = 3600 s.
In both configurations (E1 and E2), in order to compare results from both simulations and assess the resolution-dependence of the model, HR carbon concentrations are mapped to the LR discretization following Eq. (25). Moreover, in order to simplify the problem and to focus only on the resolution-dependence of the framework, coagulation and fragmentation rates, K and F 300 respectively, are set to constant values (Table 2).

Results and discussion
4.1 Linear discretization Figure 5 shows the results from the linear size discretization (E1), for both LR and HR simulations and for coagulation and fragmentation considered separately as well as simultaneously.

305
Starting with initial uniform carbon concentration distribution, coagulation leads to a reduction of C p in small size bins and an increase in larger ones for both LR and HR, resulting in a linearly increasing distribution of C p over the resolved size range (p = 1 to N , Fig. 5a,b). The largest accumulation of C p appears in the unresolved size range (p = 3 2 N ) for both LR and HR simulations.Notice that the carbon concentration in the unresolved bin of the HR experiment is divided by a factor 100 in order to fit in the y scale. However, a smaller amount of particles end up in this range in the HR simulation compared to the LR 310 Figure 5. Evolution of the carbon concentration distribution over the size range in the linear size discretization configuration (E1) as function of the arbitrary size. Our three sets of reaction simulations are represented : coagulation only (panels a, b and c), fragmentation only (panels d, e and f) and when both reactions are combined (panels g, h and i). The left column (a, d, g) represents our Low Resolution (LR) set up, the middle column (b, e, h) the High Resolution (HR) one, and the right column (c, f, i) the comparison of the HR carried back to LR (See Eq. 25 in Sect. 3.1 for details). Y-abscissas are different between our resolutions and are set to allow an easy comparison. For each panel, the initial time step t0 is in black and the final time step t1 = 24 h appears in red for LR and in green for HR. As it is the linear size discretization configuration, the LR resolved bins indexes equal the arbitrary size. Solitary points represent the size bin 3 2 N for both resolutions, as they represent the average of a larger number of size bins (LR: 5 to 8, and HR : 401 to 800) than the ones in the resolved range. Notice that the final carbon concentration in the unresolved bin of the HR experiment is divided by a factor 100 in order to fit in the y scale. As detailed in the methods (Sect. 2.3.1), both coagulation and fragmentation reactions occur in the resolved range but only fragmentation occurs in the unresolved range. For further details, refer to the legend of Fig. 4.   (Fig. 5c), which is associated with larger final concentrations C p in the resolved range for HR since the model is conservative by design. The final carbon concentration distributions are nonetheless very similar considering the very large difference in the number of bins between LR and HR.
In contrast with coagulation, fragmentation yields a reduction of C p in larger size bins to the benefit of an increase in the small ones (Fig. 5d-f). Only very small differences of C p are noticeable between LR and HR over the resolved range (Fig. 5f), 315 suggesting that fragmentation is very weakly dependent on the size range resolution when a linear discretization is used. As the unresolved range is initialized to zero, and fragmentation does not allow particle to increase in size, C 3 2 N remains zero for both simulations.
When both reactions are combined and act simultaneously, they nearly compensate each other in small size bins in both simulations (Fig. 5g,h) but to a greater extent in the HR case. For larger size bins however, fragmentation seems to operate 320 nonlinearly and more strongly than coagulation, leading to a smaller carbon concentrations related to large particles compared to the initial value. The comparison (Fig. 5i) shows a significantly greater C p in each resolved size bins for HR compared to LR, but it is the inverse for the unresolved range. This is explained by the asymmetry that exists in the mathematical formulation of coagulation and fragmentation (Eq. 9), the former being a quadratic function of concentration while the latter is a linear function. Overall, despite the fact that coagulation and fragmentation do not compensate each other, which is not a 325 prerequisite in the model, the dependence on the number bins for a given range of particle sizes is quite weak when using a linear discretization.
These results demonstrate that in a linearly size discretized configuration, the model is reasonably independent of the resolution when coagulation and fragmentation are used independently or combined (Fig. 5). They show however the importance of considering the unresolved size range in the model design, which guarantees mass conservation and unbounded growth due 330 to coagulation.

Nonlinear discretization
We now consider the results from the nonlinear size discretization model (E2), shown in Figs. 6 and 7. Recall that instead of a uniform C p distribution, the model is initialized with a distribution that is exponentially decreasing over a size range from 10 −6 m to 10 −2 m, thus spanning four orders of magnitude (Eq. (26), Table 2). In this case, results differ as a function of 335 resolution, for the reactions taken separately or combined. All LR simulations (K, F and KF) react more strongly than the HR ones.
Focusing first on the resolved range in the coagulation-only experiment, we observe a diminution of C p in the smaller size bins and an accumulation in the larger size bins in LR simulation (Fig. 6a). However in HR, this distribution pattern is limited to a small range of sizes between 10 −6 m and 10 −4 m (Fig. 6b). In addition, no significant accumulation of concentration is 340 observed in the largest size bins, and by extension neither into the unresolved range (p = 3 2 N ). This leads to a difference of C p regarding the larger size bins between LR and HR when mapped on the same grid (Fig. 6c), with the LR overestimating carbon concentration for larger size bins. Figure 7 shows the time evolution of the distributions of Fig. 6. In the LR case, the initial response (roughly 2 hours) is characterized by a fast time scale, followed by a slower response during the rest of the simulation (Fig. 7a). In contrast, in the HR case, the response is localised in the small size range and only the slower response is observed 345 (Fig. 7b). While attenuated, these biases are still clearly visible when HR is remapped on LR (Fig. 7c). Results for fragmentation only yield a similar biases between the LR and HR cases. Reactions are magnified in LR compared to HR (middle panels (d, e, f) of Figs. 6 and 7). When both reactions are combined, coagulation dominates over fragmentation to explain most of the observed changes in distribution of the LR simulation (Fig. 6g). In HR, both coagulation and fragmentation have localised effects on the smallest and largest size ranges (Fig. 6h). The comparison of HR versus LR shows a pattern that Note the logarithmic scale for the vertical axis representing particle diameter (m) and for the color scale representing carbon concentrations (mmolC m −3 ). Concentration value scales are different between our resolutions, and are set to allow an easy comparison. Dashed lines represent the bounds of each LR size bin. Dealing with size spectrum logarithmic scale, middle size bin concentration value will induce no concentration data in the lower bound of the first size bin in the LR (panels a, d and g) and for the HR mapped to LR grid (panels c, f and i).
dominates over fragmentation, which is expected for an initial concentration distribution that is highly skewed towards small particles.
These simulations demonstrate that when using a nonlinear size discretization and a nonlinear initial carbon concentration distribution, the model behaviour is significantly dependent on resolution. To attenuate this dependence, which arises from the 355 asymmetry between coagulation and fragmentation, we propose adding and tuning a penalty function that will compensate this difference as N varies.

Resolution dependency function
The residual dependence of the model to resolution in the nonlinear discretization arises mainly from the fact that the prefactors used in Eq. (22) are derived from a linear analysis, which yields an asymmetry between coagulation and fragmentation. To minimize the effect of this asymmetry, we propose multiplying both reaction terms by a resolution-dependent function f (N ) that is positive and that monotonically varying between a value to be determined at low N and 1 for N → ∞. For simplicity, we here assume an exponential function of the form with λ 5 and λ 6 positive constant parameters that were determined empirically (see Table 1).

365
Simulation results obtained with this correction factor (Figs. 8 and 9) show that both LR and HR simulations now agree much better. The LR simulation is the one that is the most impacted by this change (Fig. 8a,

Summary
We have developed a new 0D numerical model for representing coagulation and fragmentation as an interaction between three particles of arbitrary sizes. Particles are categorized in size bins that can be linearly or nonlinearly distributed along a given 375 size spectrum. In the linear configuration, E1, the total volume of suspended particulate matter (i.e. the sum of the volume of all individual particles) is also conserved. However, this is not strictly the case in the nonlinear configuration, E2, because it can happen that two particles of two different size bins can end up in the same size bin as the biggest one. By construction, the total concentration of carbon carried by particles is conserved over the resolved and unresolved range. The unique arbitrary size bin 3 2 N offers i) boundaries to avoid any exponential growth of the size range, ii) reduces falsified carbon concentration estimation 380 in the larger size bins. When absent, accumulation of particles was observed (exploratory studies of our current work which are not shown). This caveat is present in models using the sectional approach, the solution of which is to increase the number of bins, as described by Burd (2013).
Coagulation has a quadratic dependence on particle' number concentration participating to the reaction, while fragmentation has a linear dependence on the particle' number concentration. The linear configuration has a very week dependence to the size 385 spectral resolution (number of size bins N for a given size range). The nonlinear configuration has a significant dependence to resolution. This dependence can be overcome by multiplying both reaction terms with a function f (N ) such that f (N ) → 1 when N → ∞. However, further work is required to unearth what is causing the dependence, as the assumption made that reactions need to be divided by the inverse dependent resolution coefficient (Eq. 22) for both linear and nonlinear cases were partly wrong. The surmise that C p will be equally distributed between the size bins (i.e. C p = C T N ) is not respected in the 390 nonlinear case as distribution assigned is nonlinear over the size bins (Eq. (26), i.e. C p = C T ∆p ∆p ) in addition to be coupled to a nonlinear size discretization. The presented attempt to rectify this wrong assumption via the elaboration of the function (Eq. 27) for the nonlinear case allows prospects to thrive and reassure that solutions exist. Despite this, the model succeeds in representing the evolution of the size of suspended particulate matter due to the simultaneous action of coagulation and fragmentation using a low number of size bins (N = 4). This is comparable to biogeochemical models of low (e.g. Fasham et al. (1990) with seven states variables or Neumann (2000) with nine) to moderate complexity (e.g. Aumont et al. (2015) with 24 prognostic variables or Ward et al. (2012) with more than 50). But, the sensitivity of our model outcomes to many arbitrary constant parameters needs to be profoundly investigated, as some of them may be size-dependent and therefore vary along the Figure 9. Hovmöller plots of carbon concentration as a function of particle diameter over a 24 hours period for the nonlinearly discretized configuration (E2) with the application of the function which attempt to reduce the dependency to the resolution of our model (Eq. 27).
Simulations for coagulation and fragmentation considered separately (on top and middle lines respectively) and combined (on bottom line), and for LR (left column) and HR (middle column), are shown. For further details, refer to the legend of Fig. 7. size range. Such parameters are the coagulation rate and associated stickiness of particles (K, Eq. 9) and the fragmentation rate (F, Eq. 9). Other parameters values, such as the slope of the particles distribution (Ψ, Eq. 26) and parameters relying on 400 the carbon content estimation of each particles (N , Eq. 4) were chosen according to literature, determined through field or laboratory studies for specific geographic regions or ecosystem states. The total initial carbon concentration over the size range (C T ), will also require deeper investigation, as by acting jointly with other parameters it may affect reactions thresholds in the model, and the final outcomes and conclusions. Lastly, the role of coagulation and fragmentation reactions on the cell density behavior along the size range will be required (Gregory, 1997) in the parameterization of the settling velocity for each detritic 405 state variable implemented. This step will be required to incorporate the model in a 1-D environment coupled to physical fields.
Ultimately, when reliably parameterized, this model will be coupled to an upper trophic level ecological model and OGCMS that will enable addressing further questions related to the fate of particle evolution with depth.
Through our approach, a balance between a low computational cost and a proximity to particulate organic matter ecological 410 dynamics expected to be found in the ocean was fulfilled. This is a first attempt to fill the knowledge gap underlined by Boyd et al. (2019), by offering a model of particle transformations to incorporate in OGCMs. This work as well as the steps to come, will then offer new perspectives to estimate the downward carbon export in global models, as coagulation and fragmentation reactions will be characterized alongside of other known processes affecting the vertical flux of organic matter (e.g. grazing, remineralization). Thus, comparisons with previous studies will be possible to conclude on the influence to consider or not 415 these often ignored reactions on the estimation of the biological pump response to climate change.
Code availability. The current version of model and user manual are freely-available and can be downloaded from the referenced project page : http://doi.org/10.5281/zenodo.4432896 (See (Gremion and Nadeau, 2021)). The code concedes to the Creative Commons Attribution 4.0 International License (CC BY 4.0) and is written in fortran90 and requires the gfortran compiler.
Author contributions. GG and L-PN co-built the numerical framework and prepared the figures; GG ran the analysis and wrote the text, 420 GG, L-PN, CD and DD participated in discussions prior to model set up, L-PN, CD, IRS, PA and DD provided edits and reviews during manuscript preparation. DD secured funding to support the development of the model.
Competing interests. The authors declare that they have no conflict of interest.