Microstructure and composition of marine aggregates as co-determinants for vertical particulate organic carbon transfer in the global ocean

. Marine aggregates are the vector for biogenically bound carbon and nutrients from the euphotic zone to the interior of the oceans. To improve the representation of this biological carbon pump in the global biogeochemical HAMburg Ocean Carbon Cycle (HAMOCC) model, we im- 5 plemented a novel Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M 4 AGO) sinking scheme. M 4 AGO explicitly represents the size, microstructure, heterogeneous composition, density and porosity of aggregates and ties ballasting mineral and particulate organic 10 carbon (POC) ﬂuxes together. Additionally, we incorporated temperature-dependent remineralization of POC. We compare M 4 AGO with the standard HAMOCC version, where POC ﬂuxes follow a Martin curve approach with (i) linearly increasing sinking velocity with depth and (ii) temperature- 15 independent remineralization. Minerals descend separately with a constant speed. In contrast to the standard HAMOCC, M 4 AGO reproduces the latitudinal pattern of POC transfer efﬁciency, as recently constrained by Weber et al.


Introduction
Marine aggregates transfer biologically bound carbon and nutrients from the sunlit surface waters, the euphotic zone, 40 to the interior of the oceans. While uncertainty with respect to primary production estimates exists, about 4.0 to 11.2 Gt C yr −1 of biologically bound carbon is annually exported out of the euphotic zone of the global ocean (Laws et al., 2000;Najjar et al., 2007;Henson et al., 2012). The 45 net withdrawal of carbon dioxide (CO 2 ) from the ocean surface through export of carbon bound in particulate organic matter (POM) and biogenic minerals and subsequent release through microbial remineralization and dissolution during aggregate descent determine the strength of the so-called bi-50 ological carbon pump. The biological carbon pump critically depends on phytoplankton growth, the replenishment of the euphotic zone with nutrients through mixing and upwelling processes, and the efficiency of biologically bound carbon transfer from surface waters to the interior of the oceans 55 (Williams and Follows, 2011). The region and depth of carbon sequestration eventually determine the residence time of the biologically bound carbon upon recurrence at the ocean surface. Representing transport and fate of marine aggregates in Earth system models (ESMs) is therefore key in quantify- 5 ing the future evolution of biogeochemical cycles and particularly the biological carbon pump and its feedback on the Earth system under climate change (Ilyina and Friedlingstein, 2016). In the present study, we thus aim to advance the representation of marine aggregates in an ESM framework. 10 Marine aggregates are porous entities which are heterogeneously composed of POM, biogenic and inorganic minerals. The sinking velocity of marine aggregates, their microbial remineralization and zooplankton grazing govern the attenuation of vertical particulate organic carbon (POC) fluxes. 15 The sinking velocity of aggregates is primarily determined by their size. In addition, the internal microstructure, defined by the porosity and heterogeneous composition, entails high variability in excess density and thus sinking speed of aggregates (Iversen and Robert, 2015). Biogenic calcium car-20 bonate (CaCO 3 ) and opal structures, primarily formed by coccolithophores and diatoms, act as ballasting minerals in organic aggregates (De La Rocha and Passow, 2007;Armstrong et al., 2002). On the contrary, the available amount of POC, acting as glue, is suggested to limit the uptake ca- 25 pability for ballasting minerals before aggregates disintegrate Passow and De La Rocha, 2006;De La Rocha et al., 2008). Ballasting increases the POC transfer efficiency (Klaas and Archer, 2002;Balch et al., 2010;Cram et al., 2018), defined as the fraction of POC exported 30 out of the euphotic zone that reaches a particular depth, e.g., 1000 m (Francois et al., 2002). As CaCO 3 is significantly denser than opal, CaCO 3 is suggested to be a more effective ballasting material for marine aggregates (Balch et al., 2010), implying higher POC transfer efficiency in CaCO 3 - 35 production-dominated regions. Phytoplankton communities possess spatio-temporally varying patterns and prime the sinking flux ratios of detritus to ballasting minerals, i.e., the rain ratios. High CaCO 3 -to-opal ratios are found in oligotrophic regions of the mid-latitude subtropical gyres, while 40 opal is the prevalent ballasting mineral in high latitudes and upwelling-influenced equatorial regions (Balch et al., 2010). However, simple ballasting relationships on aggregates are questioned, and the prevailing plankton network is suggested as an additional driver for POC fluxes (Wilson et al., 2012;45 Henson et al., 2012;Guidi et al., 2016). For example, cell size and morphology present in the phytoplankton community are suggested as a primary determining factor for sinking velocity of marine aggregates (Laurenceau-Cornec et al., 2015;Bach et al., 2016). In turn, the attenuation of POC 50 fluxes is hypothesized to be modulated by microbial remineralization and by zooplankton grazing in oligotrophic and eutrophic regions (Guidi et al., 2009). Since temperature controls enzymatic reaction kinetics, and thus microbial remineralization of POC, slower attenuation and thus higher transfer 55 efficiency are suggested in cold high latitudes compared to warm oligotrophic regions (Marsay et al., 2015). For a long time, the aforementioned variable factors and processes, the limited understanding of aggregation and fragmentation processes that shape the aggregate size spectrum, and the sparse 60 number of data have retarded the emergence of a detailed picture of global pattern of POC fluxes attenuation and thus transfer efficiency.
However, quantification of the regionally varying POC transfer efficiency and its variability is key in understand-65 ing global biogeochemical cycles, in particular the carbon cycle (Falkowski et al., 1998). Recently, global POC fluxes have been constrained to possess high transfer efficiency in high latitudes and upwelling regions and lower efficiency in the subtropical gyres (Weber et al., 2016). The underlying 70 controls for the transfer efficiency pattern seem to exhibit a distinct latitudinal variability (Cram et al., 2018). The simplified model study of Cram et al. (2018) suggests aggregate size, ballasting of particles by CaCO 3 and opal, temperature effects on microbial aerobic and anaerobic remineralization, 75 water density, and molecular viscosity as major controls of the transfer efficiency pattern.
Processes of marine snow formation, ballasting and sinking are currently underrepresented in ESMs despite the relevance of aggregates for the transfer and sequestration of POC 80 to the deep ocean. Only a few global models explicitly incorporate aggregation of phytoplankton mechanistically (e.g., Gehlen et al., 2006;Schwinger et al., 2016) while neglecting ballasting effects or vice versa (Gehlen et al., 2006;Heinemann et al., 2019). POC sinking velocities in ESMs are typi-85 cally formulated to reproduce the Martin curve (Martin et al., 1987) or heuristically describe ballasting of POC with opal and CaCO 3 (e.g., Gehlen et al., 2006;Heinemann et al., 2019), which limits the process-based adaptation of sinking velocities under changing environmental conditions associ-90 ated with climate change.
As a first step, we develop the Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M 4 AGO) sinking scheme that explicitly represents composition, microstructure, and related properties such as poros- 95 ity and density of aggregates. We aim to consistently define marine aggregates with their in situ measurable properties in an ESM framework. We implement M 4 AGO in the global HAMburg Ocean Carbon Cycle (HAMOCC) model, which is part of the Max Planck Institute -Earth system model 100 (MPI-ESM), to explicate the emerging pattern of aggregate properties and examine their effect on sinking velocity and the global pattern of POC transfer efficiency. We particularly aim to (i) represent the POC transfer efficiency pattern of Weber et al. (2016), (ii) provide further understanding into 105 the underlying driving factors for this pattern and (iii) give insights into the impact of M 4 AGO on the global CO 2 flux pattern. We focus on the transfer efficiency pattern identified by Weber et al. (2016), as it was derived by diagnosing phosphate fluxes from World Ocean Atlas 2009 phosphate con-110 centration via inverse modeling. This approach benefits from many more observations used compared to direct flux observations (Usbeck et al., 2003;Weber et al., 2016) and can thus be regarded as, to date, more reliable than previous estimates with a partly opposing latitudinal pattern (e.g., Henson et al.,5 2012; Marsay et al., 2015).
With M 4 AGO, we represent marine aggregates at the global scale to provide a test bed for future investigations of aggregate-associated processes in ESMs, e.g., particle-size-, microstructure-and composition-dependent remineralization 10 rates.

Model description
The HAMburg Ocean Carbon Cycle (HAMOCC) model is a global biogeochemical model which features biology and resolves the carbon chemistry (Six and Maier-Reimer, 1996;Paulsen et al., 2017;Mauritsen et al., 2019). HAMOCC assumes a fixed stoichiometry for dead and living organic matter and represents the nutrients phosphate, nitrate, silicate and iron. Phytoplankton in HAMOCC, namely bulk phytoplankton and diazotrophs, can thus expe- 20 rience nutrient co-limitation. Diazotrophs assimilate gaseous dinitrogen under nitrate limitation and compete for phosphorus with bulk phytoplankton. Diazotrophs grow slower than bulk phytoplankton and have their optimal growth temperature at about 28 • C (Paulsen et al., 2017(Paulsen et al., , 2018. Zoo- 25 plankton feeds on bulk phytoplankton and releases POM, which enters the common detritus pool. During detritus formation through bulk phytoplankton or zooplankton, opal or CaCO 3 is produced depending on silicate availability. This treatment adequately depicts the spatial distribution of sili- 30 cifying and calcifying plankton communities (Heinze et al., 1999). HAMOCC represents sediment processes (Heinze et al., 1999) and is coupled to the global three-dimensional Max Planck Institute Ocean Model (MPIOM; Marsland et al., 2003;Jungclaus et al., 2013). HAMOCC is described 35 and evaluated in previous studies; for details, see, e.g., Six and Maier-Reimer (1996), , Paulsen et al. (2017), and Mauritsen et al. (2019). In the following, we therefore focus on processes in the standard version, i.e., sinking and remineralization, which we modify with the 40 M 4 AGO sinking scheme. A table with the used mathematical symbols can be found in Appendix D, Table D1.

HAMOCC standard representation of sinking fluxes and remineralization
The standard version of HAMOCC (Mauritsen et al., 2019) 45 represents sinking fluxes of POC, F POC , at depth z > z 0 , according to the concept of the Martin curve (Martin et al., 1987;Kriest and Oschlies, 2008): where F 0 is the POC flux out of the euphotic zone at ex-50 port depth z 0 . For simplicity, the export depth is globally defined as being z 0 = 100 m in HAMOCC. Above z 0 , a constant sinking speed of 3.5 m d −1 is assumed. Below z 0 , we assume a linearly increasing mass concentration-weighted mean sinking velocity with depth. The ratio between the 55 remineralization rate of POC, R POC,remin , and the vertical gradient of the sinking velocity, ∂ z w s , determines the POC flux slope, β = R POC,remin /∂ z w s (Kriest and Oschlies, 2008).
In the standard version of HAMOCC, remineralization of POC is temperature-independent and comprehends oxygen-60 concentration-dependent aerobic remineralization as well as sulfate reduction and denitrification under sub-anoxic and anoxic conditions. The sinking tracers, opal and CaCO 3 , are treated separately from POC and sink with their own, constant sinking 65 velocity. Aeolian dust is, apart from the release of bioavailable iron in surface waters, inert and sinks slowly through the water column. The opal dissolution rate in the standard model is linearly temperature-dependent. HAMOCC accounts for dissolution in carbonate ion undersaturated con-70 ditions below the dynamically emerging lysocline. In the following, we refer to this version of HAMOCC as "standard".

The novel M 4 AGO sinking scheme in HAMOCC
Natural waters exhibit a size spectrum of aggregates whose diameter, d, composition and microstructure determine their 75 terminal sinking velocity. In the M 4 AGO approach, we explicitly represent microstructure and heterogeneous composition of aggregates. For the aggregate size spectrum, we limit the representation to a variable power law number distribution, n(d), with slope b and power law factor a (fol-80 lowing e.g., Kriest and Evans, 1999;Gehlen et al., 2006;Schwinger et al., 2016): This way, we avoid the computational costs of size-classbased model approaches (Jackson, 1990;Stemmann et al., 85 2004;Sherwood et al., 2018). The local concentration-weighted mean sinking velocity, w s , in M 4 AGO is eventually computed from the number distribution (Eq. 2) that is truncated at the minimum and maximum aggregate sizes, d min and d max , respectively, and 90 expressions for the aggregate mass, m(d), and the sinking velocity of aggregates, w s (d), of a particular diameter, d. Integration over the aggregate size spectrum yields w s We only implicitly account for aggregation and fragmen-95 tation and explicitly represent the temporally and spatially www.biogeosciences.net/17/1/2020/ Biogeosciences, 17, 1-39, 2020 variable heterogeneity and microstructure of aggregates and their effect on the mean sinking velocity. We refrain from representing the potential heterogeneity of aggregate composition within the local size spectrum (see, e.g., Jackson, 1998). Consequently, and in contrast to the standard configuration, the settling tracers in HAMOCC, opal, CaCO 3 , detritus and dust sink in M 4 AGO at the same mean sinking velocity of aggregates (Eq. 3). In contrast to Gehlen et al. (2006), Schwinger et al. (2016)  Marine aggregates are porous ) and feature a self-similar microstructure which can be described via a fractal dimension d f (Logan and Wilkinson, 1990;Kranenburg, 1994). d f = 1 would depict a chain of aggregate constituents, where the length equals the aggregate 20 diameter, and d f = 3 describes a solid sphere. Consequently, the mass of an aggregate, m(d), grows disproportionately to the aggregates volume and can be expressed as where m f is a mass factor for the smallest entity. Thus, the 25 density of an aggregate ρ f decreases with increasing diameter. Aggregates consist of, for example, phytoplankton cells or coccolithophore shells , which we consider to be spherical primary particles. Primary particles exhibit their own density ρ p and diameter d p . Tak- 30 ing the fractal scaling of aggregate mass into account, the excess density of an aggregate ρ f = ρ f − ρ with respect to surrounding fluid density ρ can be expressed as (Kranenburg, 1994)

35
Furthermore, the aggregate porosity, φ, is defined as and, hence, both excess density and porosity are regulated by the fractal dimension and primary-particle size. The ρ f can be introduced to the well-known Stokes (1851) terminal 40 sinking velocity, w s : For small particle Reynolds numbers, Re p = w s d/ν < 0.1, the drag coefficient is c D = 24/Re p and the sinking velocity, w s , becomes where µ and ν are molecular dynamic and kinematic viscosity (Matthäus, 1972), respectively, and g is the gravitational acceleration constant. However, this approach assumes homogeneous, mono-sized primary particles, while it displays 50 the potential importance of primary-particle size and density as well as aggregate microstructure for sinking velocity.
To better represent aggregates in natural systems, the heterogeneity of primary particles was thus far considered either for size or density (Jackson, 1998;Maggi, 2009;Khelifa 55 and Hill, 2006). With M 4 AGO, we represent aggregates composed of poly-dense, poly-sized primary particles under the assumption of a singled value fractal dimension throughout the aggregate size spectrum. This allows for representing heterogeneous primary particles such as diatom frustules, coc-60 coliths, dust particles and detritus as principal components of marine aggregates. Bushell and Amal (1998) derived a representation of the mean primary-particle size, for an aggregate that is composed of n p = i n i mono-dense spherical primary particles of different diameters d p,i . Polysized formed aggregates disobey the traditional mass fractal relationship, but the fractal nature continues to emerge in a power law scaling for the mass present in a radial shell from 70 an occupied point in the aggregate (Bushell and Amal, 1998). The approach of Bushell and Amal (1998) conserves the size of the aggregate and the encapsulated solid volume of the primary particles, and thus the porosity of the aggregate is unimpaired, while the calculation does not presume an equal 75 number of mean, n, and individual primary particles, i n i (hence, n · d p 3 = i n i d 3 p,i with n = i n i for poly-sized primary particles), which is negligible in the following, as we do not consider n any further.
Under the assumption that aggregates feature the same 80 composition and hence same heterogeneity in a size spectrum, the aggregate-composing primary-particle types are always the same for any aggregate of diameter d in a unit volume and thus n i /n p = constant. This further implies that the ratio K i , between the total number of primary particles of 85 one particle type, n i,tot , and the total number of primary particles, n i,tot , in a unit volume is equal to the ratio found in an individual aggregate: Rewriting Eq. (10), n i = K i i n i , and inserting it into Eq. (9)+ gives where the factors K i can be expressed in HAMOCC via the concentration of each aggregate-forming tracer C i . Namely, 5 we consider the HAMOCC tracers detritus, opal, calcite and dust in taking part in the formation of heterogeneously composed aggregates. Calculating the number of primary particles from the tracer concentration requires the molecular concentration to mass factor, R i , the tracer-related primary-10 particle diameter, volume V p,i = 1 6 π d 3 p,i , and density ρ p,i : The advantage of Eq. (11) is that it allows us to determine the mean primary-particle diameter in HAMOCC while solid volume and density of primary particles are conserved. En- 15 suring mass conservation, we introduce the volume-weighted primary-particle mean density and, hence, multiplication by the volume of the mean primary particle then yields the mass of a mean primary particle 20 (see also Fig. 1): Substituting Eq. (14) into Eq. (4), we derive the mass factor for heterogeneous aggregates m f = 1 6 π d p 3−d f ρ p . The derivation of the mean primary-particle diameter (Eq. 11), 25 density (Eq. 13) and mass (Eq. 14) allows for applying common fractal laws for the calculation of aggregate mass, density and thus sinking velocity. Hence, w s (d, ρ p , d p , . . .) can be expressed as w s (d, ρ p , d p , . . .). For a single type of primary particle, all underlying equations reduce to the tra-30 ditional fractal-scaling relationship (Logan and Wilkinson, 1990;Kranenburg, 1994).

Mean sinking velocity of marine aggregates
In the preceding section, we derived a formulation for the mean primary-particle size (Eq. 11), which we apply as a 35 lower integration bound in Eq. (3), and, hence, d min = d p (following Kriest and Evans, 1999). The maximum aggregate diameter of the size spectrum, d max , is limited by fragmentation of particles. Several mechanisms can cause fragmentation of aggregates. Flow-induced turbulent shear has 40 been suggested as the dominant process in the upper ocean, where turbulent shear reaches typical values of the order of Figure 1. Underlying assumptions for the representation of aggregates composed of poly-dense, poly-sized primary particles. Primary particles, like dust particles, coccoliths and diatom frustules (a), are assumed to be spherical and exhibit their characteristic density (b). Once aggregated, we assume the diameter of the aggregate to be constant and the total volume and mass of primary particles to be preserved (c). V p,i is the primary-particle volume, ρ p,i is the primary-particle density and d p,i is the primary-particle diameter of primary-particle type i. m(d) is the mass of an aggregate of diameter d. d p and ρ p represent mean primary-particle diameter and density, respectively. 1 s −1 (Jackson, 1990). By contrast,  showed that marine aggregates often withstand oceanic turbulence conditions and suggested biological processes as a 45 mediating factor for shaping the size distribution. Zooplankton also generates turbulent shear that is strong enough to rupture aggregates (Dilling and Alldredge, 2000). Hill (1998) proposed an alternative control on aggregate size, namely the sinking of aggregates that produces shear of the same order 50 of magnitude as ambient turbulent shear in the ocean (Bagster and Tomi, 1974;Adler, 1979;. Sinking could thus cause fragmentation in deeper regions of the ocean, where turbulent shear is small (O(0.01 s −1 ); McCave, 1984;Waterhouse et al., 2014). Since modeling of 55 particle-reactive thorium suggests continued fragmentation during particle descent in the deep ocean (Lam and Marchal, 2015), we adopt the hypothesis of sinking-induced fragmentation and limit the size distribution based on the particle Reynolds number, Re p ,: This drag representation leads to smaller settling velocities for large aggregates than the classical Stokes drag (c D = 24/Re p ). Hence, aggregates can grow larger until they reach the globally fixed critical Re p for fragmentation, Re crit , which leads to a more realistic representation of the size range of aggregates. We approximate the White drag representation to be (Jiang and Logan, 1991) to avoid iteration and to allow for an analytical solution of Eq. (3). Applying the parameter values of a j =1 = 24.00, b j =1 = 1, for Re p ≤ 0.1; a j =2 = 29.03, b j =2 = 0.871, for 0.1 < Re p ≤ 10; and 10 a j =3 = 14.15, b j =3 = 0.547, for 10 < Re p ≤ 100 introduces maximum errors of less than 10 % compared to Eq. (16) for Re p < 100 (Jiang and Logan, 1991).
By introducing Eq. (17) in Eq. (7), and applying Eq. (5) using mean primary-particle properties, the approximation for 15 the sinking velocity becomes Consequently, the concentration-weighted mean sinking velocity (Eq. 3) can then be expressed as

25
where d max = d j (Re crit ) is the maximum diameter of aggregates, and by applying a j =0 = b j =0 = 1, the lower integration boundary equals the mean primary-particle diameter.

The particle distribution slope, b
Observed aggregate size spectra in the ocean exhibit a spatio- 30 temporally dependent slope ranging between approximately 3.2 and 5.4 (DeVries et al., 2014) or even lower (≈ 2; Guidi et al., 2009). A smaller slope parameter, b, translates to more large aggregates relative to a larger b and enhances mean sinking velocity. The evolution of the particle-size spectra 35 underlies the interacting processes of growth and decay of phytoplankton, aggregation, fragmentation and sinking of aggregates. Instead of modeling the processes of aggregation and fragmentation explicitly or prescribing b, we assume a dynamic steady state between aggregation and fragmentation 40 to describe the slope of the number distribution. According to dimensional analysis, the slope of the number distribution in dynamic steady state depends on the fractal dimension of aggregates and the process of aggregation, aggregation due to shear, differential sinking and Brownian motion (Jiang and 45 Logan, 1991). Aggregation due to Brownian motion is only relevant for particles smaller ≈ 1 µm (McCave, 1984), which we neglect here. We further assume that aggregation in the majority of the global ocean is dominated by differential settling and express the particle distribution slope, b, as (Jiang 50 and Logan, 1991) where b J is a fixed parameter for the sinking velocity dependency on the particle Reynolds number that we fix for simplicity to b J = b j =2 . The assumption of differential-settling-55 dominated aggregation is likely violated in the euphotic zone, where shear aggregation is probably more relevant and steady-state assumption is questionable, which we will address in the discussion (Sect. 3.10).

stickiness and fractal dimension
Adhesion properties of particles affect the fractal structure of aggregates and the collision efficiency ("stickiness") of particles (Meakin, 1988;Liu et al., 1990). Theoretical studies show that the stronger the surface adhesive forces are, 65 the higher the stickiness of particles and the smaller the intrusion of particles and particle clusters into each other are (Liu et al., 1990). As a result, this leads to a looser structure, which translates to a small fractal dimension. Stickiness of phytoplankton is species-specific (Hansen and Kiørboe, 70 1997) and depends on the growth phase (Simon et al., 2014). Furthermore, phytoplankton releases extracellular polymeric substances (EPS; Decho, 1990) such as transparent exopolymer particles (TEPs), which are suggested to be aggregationpriming, sticky materials (Alldredge et al., 1993;Azam and 75 Malfatti, 2007;Thornton, 2002;Passow, 2002;Engel et al., 2004). The resulting fractal dimension is typically determined as one value across all aggregate sizes. We thus assign a single fractal dimension to an aggregate population and depict the linkage between stickiness and fractal dimension 80 in a qualitative manner. We attribute a stickiness value, α i , to each sinking tracer in HAMOCC and calculate the mean stickiness for the aggregates that is then mapped to a fractal dimension. Since adhesion, and thus stickiness, is a surface property, we calculate the mean stickiness, weighted by the primary-particle surfaces A i ∝ d 2 p,i . We map the mean stickiness to a range between 0 and 1: where α min = min(α i ) and α max = max(α i ). 15 , where d f,min and d f,max are the parameterized minimum and maximum fractal dimension of aggregates. Modeled stickier aggregates thus exhibit lower fractal dimensions than non-sticky particles, which is in qualitative agreement with previous studies (Meakin, 1988;Liu 20 et al., 1990;Block et al., 1991;Nicolás-Carlock et al., 2016).

Diatoms as a special case of primary particles
Diatoms are silicifying phytoplankton that possesses a hollow opal skeleton, the diatom frustule, and are thus different from a homogeneous, solid primary particle like a coc- 25 colith. Diatoms feature a wide range of sizes, from about a few microns to millimeters (Armbrust, 2009). Since sinking velocities of aggregates are proportional to their diameter, primary-particle density and size (Eq. 8), aggregateincorporated large diatom shells likely enhance the sinking 30 velocity of particles. Indeed, un-remineralized, intact diatom frustules were even found in deeper regions of the ocean (Assmy et al., 2013), which is in agreement with previously found high sinking speeds of large diatom aggregates (Alldredge and Gotschalk, 1988). We therefore explicitly account 35 for diatom shells by treating them as hollow opal spheres, filled (i) with detritus and (ii) increasing water content with ongoing remineralization while sinking (see Fig. 2). The opal volume of a modeled diatom is

40
where d p,fustule is the diameter of the diatom, whose opal shell thickness l is expressed in terms of the fixed opal-tophosphorus formation ratio. The number of diatom frustules per unit volume, can therefore be deduced from the present opal concentration, [opal], and the opal mole-to-weight factor R opal according to Eq. (12). We assume that the modeled detritus pool can be split into a free external, non-diatom and a diatom frustule-related, void-filling detritus part. We further assume 50 that the external pool is remineralized before the intracellular pool of volume V POM and thus neglect cell lysis observed prior to aggregation (Armbrecht et al., 2014) and rather assume mineral protection of detritus (Hedges et al., 2001). If more detritus is remineralized than the frustules void would 55 hold, it is replaced with the respective volume of water V aq of density ρ. The frustule density is thus During growth and decay, diatoms excrete TEPs which are positively buoyant and possess a density of about ρ TEP = 700 60 to 840 kg m −3 (Azetsu-Scott and Mari et al., 2017). TEPs are suggested to play a prominent role in aggregation processes, as they are probably sticky and thus enhance aggregation (Dam and Drapeau, 1995;Passow, 2002). In HAMOCC, phytoplankton excretion of TEPs is not re-65 solved explicitly. We therefore treat TEPs virtually and assume a linear dependency of diatom stickiness and density on the freshness of detritus, defined as the mass ratio between the actual amount of detritus, m e = n frustule V POM ρ det , and the potential mass of detritus linked to diatom frustules, 70 m potential = n frustule (V POM +V aq ) ρ det . An additional underlying assumption is that TEPs are remineralized with the same rates as normal detritus. Hence, we define the stickiness of diatoms as for m potential > 0, where α TEP and α opal are the stickiness of TEPs and pure opal, respectively. To account for the additional buoyancy through TEPs (Jokulsdottir and Archer, 2016;Mari et al., 2017), here we simplify and assume that the frustule density is lowered by TEPs in dependency on the 80 8 J. Maerz et al.: Microcomposition of marine aggregates as co-determinant for global POC fluxes freshness of detritus. Eventually, the diatom density, ρ diatom , becomes TEPs thus have a 2-fold effect on aggregates in our model: (i) TEPs increase stickiness and loosen the aggregate structure, thus decreasing the fractal dimension of aggregates, and (ii) TEPs decrease the fresh diatom frustules' density and thus add buoyancy without violating tracer mass conservation (see also model discussion in Sect. 3.10).

Temperature-dependent opal dissolution and POC
10 remineralization Marine aggregates tie heterogeneous components together that are disparately remineralized or dissolved. By contrast, in the standard model, detritus, opal and CaCO 3 were sinking separately from each other, and the global remineraliza-15 tion and dissolution rates are tuned independently because the processes are artificially decoupled. In M 4 AGO, remineralization of detritus and dissolution of, in particular, opal is tightly linked through the same sinking velocity which let us re-evaluate and revise the formulations for opal dissolution 20 and remineralization. Opal dissolution is temperature-dependent (Ragueneau et al., 2000 and is microbially mediated (Bidle et al., 2002). Intact diatom frustules are protected from dissolution by an organic matrix (Lewin, 1961). Once the organic 25 protection surrounding the silicate frustule becomes utilized by temperature-dependent microbes, they initiate and mediate the dissolution of opal (Bidle et al., 2002). Hence, opal dissolution follows a sequential process: (i) an initial temperature-dependent remineralization of the organic coat- 30 ing of the silicate frustule and (ii) the microbially mediated dissolution of opal with a temperature dependency of Q 10 ≈ 2.3 (Bidle et al., 2002). We here focus on the temperaturedependent microbially mediated dissolution and introduce a Q 10 temperature-dependent opal dissolution: where r opal is the opal dissolution rate at the reference water temperature T ref,opal and T is the ambient water temperature. In the standard version, we remain with the former linearly temperature-dependent opal dissolu- Ragueneau et al., 2000;Segschneider and Bendtsen, 2013). Analogously to opal, we incorporate a temperaturedependent Q 10 factor to aerobic POC remineralization (Dell et al., 2011;Mislan et al., 2014) We keep the anaerobic remineralization temperatureindependent, since we do not expect temperature shifts in ocean depths, where oxygen minimum zones appear in HAMOCC. In the standard run, the remineralization rates are temperature-independent (Q 10,POC = 1).  -40;Simmons and Gibson, 2000;Röske, 2005). The mean annual cycle of wind stress, heat and freshwater fluxes is 70 resolved on a daily basis. The continental freshwater runoff is provided by means of a runoff model (Röske, 2005). The loss of POM, opal and CaCO 3 due to sedimentation and subsequent burial was accounted for through homogeneously applied weathering rates which were adjusted for 75 the standard run (and the M 4 AGO run): globally, we add ≈ 99.6 (101.5) G mol P yr −1 as dissolved organic phosphorus and ≈ 3.2 (2.3) T mol Si yr −1 . To compensate for the loss of CaCO 3 , we add ≈ 17.2 (26.5) T mol C yr −1 to surface dissolved inorganic carbon (DIC) and a corresponding 80 amount to surface total alkalinity, A T , as in DIC : 2A T . We start the M 4 AGO run from the standard run at steady state and spin it up until steady state is reached in surface and mesopelagic waters, which translates to 1700 model years.
Through the long overturning times of the global ocean, we 85 still see drifts of nutrient concentrations in deep, old North Pacific waters at this state (i.e., on average ∼ 7.1 µmol P m −3 per century below 2000 m, which amounts to a centennial change of about 0.25 %). We neglect this drift, as we focus on the aggregate properties and their effects on POC fluxes 90 throughout the euphotic and mesopelagic zone.

Parameters of the M 4 AGO scheme
The M 4 AGO sinking scheme introduces a set of new parameters, in particular the primary-particle characteristics, which require constraining and tuning (summarized in Table 1). 95 We applied HAMOCC's standard sediment densities for opal, CaCO 3 and dust to the densities of primary parti-  Primary particles featuring size, density and stickiness are detritus, diatom frustules, coccoliths (CaCO 3 ) and dust minerals. The number of diatom frustules (Eq. 26), related diatom density (Eq. 29) and stickiness (Eq. 28) are estimated from opal and detritus concentration. Diatoms are then considered to be primary particles which feature particular characteristics. The remaining detritus is considered to be detritus primary particles. The calculation of the fractal dimension (Eq. 24) and the calculations of mean primary-particle size (Eq. 11) and density (Eq. 13) are carried out. The fractal dimension determines the number distribution slope (Eq. 21).
The minimum, d min = d p , and maximum, d max , aggregate diameter (Eq. 19) are estimated.
The mean sinking velocity (Eq. 20), with which the tracers sink, can eventually be determined. cles, namely ρ opal = 2200 kg m −3 , ρ calc = 2600 kg m −3 and ρ dust = 2600 kg m −3 . For suspended detritus, we chose ρ det = 1100 kg m −3 (Fettweis, 2008). As density of TEPs, we applied ρ TEP = 800 kg m −3 , which is within the measured range of 700 to 840 kg m −3 (Azetsu-Scott and Passow, While the density of primary particles is comparably well constrained, the adhesion forces of primary particles, namely related stickiness and fractal dimension of aggregates, are less studied and are weakly constrained. There is still no 10 standardized way to investigate stickiness, fractal dimension and their interdependence for aggregates in natural waters. Stickiness is experimentally defined as the interparticle attachment rate divided by the interparticle collision rate. Uncertainties in either of the two rates, e.g., due to ignoring the fractal structure of aggregates, aggregate permeability, etc., directly affect the calculated stickiness (Filella, 2007). In addition, stickiness is phytoplankton species-specific (Hansen and Kiørboe, 1997) and depends on the growth phase (Simon et al., 2014). The methodological limitations, the het-20 erogeneity of marine aggregate constituents and their variable formation process lead to a wide spread of indirectly inferred values for stickiness and fractal dimension (see, e.g., Filella, 2007, for a broader overview). Diatom aggregates seem to feature a wide spread of fractal dimensions ranging 25 from d f ≈ 1.26 to d f ≈ 2.46 (Alldredge and Gotschalk, 1988;Guidi et al., 2008) and α of about 0.03 to 0.88 (Kiørboe et al., 1990;Dam and Drapeau, 1995;Alldredge and McGillivary, 1991). Indirectly inferred fractal dimensions for reworked aggregates exhibit values of d f ≈ 2.26 to d f ≈ 2.46 30 (Guidi et al., 2008), and mineral-dominated aggregates also feature high fractal dimensions of about d f ≈ 2 ( Winterwerp, 1998) to d f ≈ 2.6 (Kranenburg, 1999) and low stickiness of α ≈ O(10 −2 ) (Tambo and Hozumi, 1979). Since stickiness of in situ primary particles is seldom measured, 35 we choose it to our best knowledge and order the stickiness for modeled primary particles according to the mean stickiness of observed aggregate types: α dust < α opal ≤ α CaCO 3 < α det < α TEP (see also Table 1). Hence, detritus and TEP-rich aggregates are modeled with a loose structure and low frac-40 tal dimension, while degraded, mineral-rich aggregates are more compact and thus feature a higher fractal dimension (see Eq. 24) which is in congruence with the present conceptual understanding of aggregates becoming compacted during their descent (Mari et al., 2017). For the minimum and 45 maximum fractal dimension of marine aggregates, we chose conservative bounds of 1.6 (Logan and Alldredge, 1989;Li and Logan, 1995;Alldredge, 1998) and 2.4 (in the range of d f ≈ 2.26 to 2.46 for reworked aggregates; Guidi et al., 2008), which is well within the observed range of 1.26 (Lo-50 gan and Wilkinson, 1990) to 2.6 (Kranenburg, 1999) for marine particles.
For the primary-particle sizes, we conceptually assume that the tracer characteristics of opal and CaCO 3 are primarily related to phytoplankton mineral structures such 55 as diatom silicate frustules and the coccoliths of coccolithophores. This implies that modeled zooplankton egests biogenic mineral structures of algae, while their own larger mineral body structures play, in numbers, a minor role in biogenic mineral fluxes (as described by, for example, Berelson 60 et al., 2007;Ziveri et al., 2007;Fischer and Karakas, 2009;Fischer et al., 2016). Coccolith diameters range from about 1.5 to 15.5 µm (Young and Ziveri, 2000;Henderiks and Pagani, 2008;Henderiks, 2008). The globally ubiquitous coccolithophore Emiliania Huxleyi (Read et al., 2013) exhibits 65 coccoliths of about 3 to 4 µm in diameter (Young and Ziveri, 2000). We set d p,calc = 3 µm, which is thus at the lower bound of the observed range, to account for the volumetric density effect of non-spherical plate-like coccoliths. Diatom frustules feature sizes of a few micrometers to millimeters 70 (Armbrust, 2009), with a dominant size range of about 12 to 58 µm (equivalent spherical diameter of body volumes 10 3 -10 5 µm 3 ; Litchman et al., 2009). We define the diatom frustule size as d p,frustule = 20 µm. The primary-particle size of detritus is weakly constrained and likely ranges from sub-75 micrometers of microgels and bacteria to millimeter scales of zooplankton body structures (Verdugo et al., 2004). We here chose d p,det = 4 µm. Aeolian dust particles feature a typ-ical size ranging from submicron to about 20 µm in size, with a mass median diameter of about 1.5 to 3 µm (Maher et al., 2010;Mahowald et al., 2014). In summary, we assumed the following order of primary-particle sizes for the tracers: d p,dust < d p,calc ≤ d p,det < d p,frustule .
The size distribution-limiting maximum aggregate diameter, d max , is variable in the model domain and depends on the critical particle Reynolds number Re crit . Re crit and its potential dependency on aggregate properties are weakly constrained, which lets us fix the value globally to Re crit = 20, 10 a conservative value compared to the measured maximum particle Reynolds number of up to 32 by Alldredge and Gotschalk (1988).
Opal dissolution and detritus remineralization are Q 10 temperature-dependent in M 4 AGO. Typically, the Q 10 fac- 15 tor for biological processes is in the range between 2 and 3. We here chose Q 10,POC = 2.1 (similar to Mislan et al., 2014, who applied Q 10,POC = 2.0). For opal, we tuned Q 10,opal = 2.6 compared to 2.3 suggested by Bidle et al. (2002). 20 The newly parameterized processes of sinking and remineralization directly affect the transfer efficiency and thus the climatological nutrient fields. This close connectedness hampers the clear distinction between data employed for model tuning or for model evaluation, when comparing the 25 model results to literature values for transfer efficiency (Weber et al., 2016) and World Ocean Atlas data (Garcia et al., 2014a, b). The transfer efficiency, in combination with the general circulation pattern, affects the nutrient climatology in the long term. In turn, sinking velocity, remineralization 30 and dissolution define the transfer efficiency on timescales of days to months. A direct comparison of modeled to observed sinking velocities and fluxes is challenging, as scale dissimilarities introduce uncertainty for comparisons between models and observations (Bisson et al., 2018). Furthermore, sedi-35 ment trap data for POC and mineral fluxes exhibit high uncertainties which complicate model comparisons and even make different parameterizations for vertical fluxes undistinguishable (Cael and Bisson, 2018). As a consequence, we perform a general evaluation of our model results. 40 Long simulations with high computational costs to reach steady state are required in the process of model tuning, which prevents intensive parameter variations. We performed parameter variations aiming at a quantitative agreement with the transfer efficiency of Weber et al. (2016). Since the ad-45 justment of the sinking velocity versus the remineralization and dissolution rates, and thus the transfer efficiency, occurs within a few years, this strategy was useful for selecting for promising parameter sets. With respect to the primaryparticle characteristics, we kept the stickiness values, once 50 chosen to our best knowledge, untouched and minimally varied the primary-particle sizes within the range of literature values. We primarily focused on tuning the remineraliza-tion and dissolution rates of POC and opal, respectively. We choose this strategy since reliable remineralization and disso-55 lution rate measurements are available to evaluate the tuned rates (see ranges in Table 1). We aimed at keeping global mean values of primary production, export production of POC, opal and CaCO 3 , and their fluxes to sediment within estimated literature ranges. This let us minimally vary the 60 fraction of maximum CaCO 3 production, c max , in M 4 AGO compared to the standard run.

Model tuning and evaluation
We here compare and evaluate the M 4 AGO run with respect to (i) where possible, the standard run, (ii) the regional transfer efficiency, as derived by Weber et al. (2016), 65 (iii) World Ocean Atlas data from the World Ocean Database (Boyer et al., 2013;Garcia et al., 2014b), and, independently, (iv) sediment trap-sampled POC and biogenic mineral fluxes compiled by Mouw et al. (2016a, b). The Mouw et al. (2016a, b) data were time-and depth-weighted to re-70 ceive monthly climatological values for the respective grid boxes in MPIOM, where the sediment trap records were taken. Model results are presented as yearly mean of the last simulated year, unless stated otherwise.

75
In the following, we evaluate the global net primary production, the export of POC to the mesopelagic zone and the associated pattern of biogenic mineral fluxes (Sect. 3.1). In M 4 AGO, the pattern of POC and associated minerals determines the aggregate properties, which we explicate in 80 Sect. 3.2. In Sect. 3.3, we present the global pattern of transfer efficiency. In Sects. 3.4 and 3.5, we examine the contributions of remineralization rates, sinking velocity and aggregate properties to the transfer efficiency pattern. Thereafter, we evaluate the modeled rain ratios (Sect. 3.6) and the bio-85 geochemical tracer distributions (Sect. 3.7). In Sect. 3.8, we discuss the consequence of the transfer efficiency pattern on regional CO 2 fluxes. Subsequently, we examine the sensitivity of the transfer efficiency to selected model parameters (Sect. 3.9) and conclude with a critical review of underlying 90 assumptions of M 4 AGO (Sect. 3.10).

Spatial distribution of POC export fluxes and associated biogenic minerals
The global pattern of the depth-integrated primary production is dominated by global circulation and thus nutrient 95 transport. The global pattern of annual mean integrated primary production therefore remains similar between the standard and the M 4 AGO run (Fig. 4). Globally integrated, the annual net primary production is ≈ 55.3 Gt C yr −1 in M 4 AGO compared to ≈ 44.7 Gt C yr −1 in the standard run. 100 The ratio between carbon export out of the euphotic zone at depth z 0 = 100 m and the net primary production, NPP, provides an estimate of how efficient the export is with respect to the net primary production. Globally, about 5.56 5 and 6.28 Gt C yr −1 is exported out of the euphotic zone in M 4 AGO and the standard run, respectively. In the standard run, the subtropical gyres exhibit p ratios of more than 0.2, and the high latitudes feature lower export efficiency (Fig. 4c). In the M 4 AGO run, the equatorial Pacific exhibits the lowest export efficiencies, the subtropical gyres feature maximum values of about 0.14-0.16 and the Arctic region maximum value is about 0.20 (Fig. 4d). The M 4 AGO run thus possesses smaller latitudinal variability in the p ratio compared to the standard run. 15 In comparison to previous estimates on global primary production, both model runs are well within the range of 30 to 70 Gt C yr −1 and show similar pattern of NPP (Carr et al., 2006). The higher NPP in M 4 AGO is due to the enhanced remineralization rates in surface waters, which also lead to 20 the lower export efficiencies in the equatorial and subtropical regions. Estimates of the export efficiency from satellite data, in situ observations or models lead to partly contrasting patterns (e.g., Lutz et al., 2002;DeVries and Weber, 2017;Henson et al., 2011Henson et al., , 2012Buesseler, 1998;Neuer et al., 2002;25 Siegel et al., 2014;Cram et al., 2018). In contrast to our two model runs, the highest p ratios are suggested to be found in the North Pacific and Antarctic Ocean (e.g., Dunne et al., 2005;Henson et al., 2011;DeVries and Weber, 2017). In the tropical and subtropical regions, the M 4 AGO run reduces the bias with respect to previously found low export efficiencies of about 1 % to 10 % (Buesseler, 1998;Neuer et al., 2002).
The exported POC is accompanied by biogenic minerals and dust. The export flux ratios between opal and detritus ex- hibit a clear latitudinal pattern, with high values in the high 35 latitudes and upwelling regions compared to the subtropical gyres (Fig. 5a, b). Both model simulations show a similar pattern of opal-to-detritus flux ratios. Higher CaCO 3 -todetritus flux ratios are confined to equatorial and subtropical regions where silicate depletion favors calcification (Fig. 5c, 40  d). M 4 AGO exhibits a higher CaCO 3 -to-detritus-mass flux ratio in the western tropical and subtropical Pacific than the standard run. The higher remineralization in the surface waters in this region reduces the amount of detritus that can coalesce with CaCO 3 . The spatial distributions of sinking trac-45 ers and their ratios prime the marine aggregate properties in M 4 AGO.

Spatial distribution of marine aggregate properties
The spatial patterns of detritus and mineral fluxes are reflected in the distribution of aggregate properties (cf. pat-50 tern in Figs. 5b and d to 6a-f). The primed characteristics further evolve while aggregates descend through the water column and become remineralized. Hence, the information of the tracer distribution in the euphotic zone propagates into the mesopelagic zone. The mean primary-particle 55 density, ρ p , ranges at export depth (100 m) from about 1100 kg m −3 in diatom-dominated regions to 1850 kg m −3 in the western Pacific, where CaCO 3 -to-detritus export ratios are high (Figs. 5b and d and 6a). In the Arctic, some regions harbor ρ p of maximum 2600 kg m −3 where aggre-60 gates in our model are dominated by dust particles (Fig. 6a). Particularly in regions of high CaCO 3 -to-detritus export ratios, ρ p increases with depth ( Figs. 6a and g and 7a). Accordingly, the volume-weighted mean excess aggregate density, ρ f V , exhibits the same pattern as ρ p , while it 65 ranges from about 2 kg m −3 in diatom-dominated regions to 35 kg m −3 in calcifier-dominated regions (Fig. 6d). Note that we chose ρ f V to account for the increasing poros-ity with size (Eq. 6), thus decreasing aggregate excess density with size (Eq. 5). As a mean value, ρ f V thus underestimates excess density for small aggregates, while it overestimates excess density for large aggregates. Generally both ρ p and ρ f V tend to increase with depth (Fig. 7a,  d). In the Pacific, CaCO 3 as ballasting mineral becomes dissolved below the lysocline, and modeled ρ p decreases again in the deep ocean (Fig. 7a, d). The mean primaryparticle size, d p , ranges between the attributed minimum and maximum primary-particle size of tracers, 2 to 20 µm, 10 respectively. Mean primary-particle size shows an opposing pattern to ρ p (Figs. 6a, b, g and h and 7a and b), since regions are either dominated by small, dense coccoliths or large, less dense diatom frustules. The global pattern of d p , primed through export fluxes of detritus and minerals, varies 15 only little throughout the water column (Fig. 7b), since we assumed invariance of primary-particle size to remineralization, dissolution and other processes. Our simulated d max is largest in the surface waters of the Southern Ocean and upwelling regions, where TEP-rich ag-20 gregates prevail. In calcifier-dominated regions, the maximum aggregate size is small. Generally, d max tends to decrease with depth (Fig. 7e).
The microstructure of marine aggregates, modeled as fractal dimension, d f , shows a pronounced spatial distribution. At 25 100 m depth, d f ranges from about 1.7 in upwelling regions and the Southern Ocean to about 2.2 in the western equatorial Pacific and features maximum values of 2.38 in Arctic regions (Fig. 6c). With increasing depth in the mesopelagic zone, aggregates tend to experience a rapid compaction as 30 d f increases (Fig. 7c). Ongoing POM remineralization during aggregates' descent shifts the aggregate composition towards mineral components, which feature lower stickiness in our model. Thus the fractal dimension increases, which mimics compaction of aggregates. The global pattern tends 35 to homogenize with depth at about 1000 m, where modeled aggregates feature d f values of about 2.2 (Fig. 6i). Exceptions are upwelling regions, where d f remains low, since detritus is slowly remineralized anaerobically in the associated oxygen minimum zones (OMZs). Below 1000 m, d f only increases 40 slowly with depth (Fig. 7c).
Particle properties and molecular dynamic viscosity determine the concentration-weighted mean sinking velocity of aggregates, w s (Fig. 6f, l). For w s , M 4 AGO considers particle sizes ranging from few micrometers to millimeters  (Fig. 6f, l). Upwelling-influenced surface waters tend to show smaller w s than the western equatorial Pacific. Generally, w s appears to increase rapidly with depth 55 within the mesopelagic zone (Fig. 7f). At about 1000 m, the latitudinal pattern changes and the high latitudes, in the Southern Ocean at around 45 • S, exhibit the highest w s of about 65 m d −1 . Along the Equator, outside the OMZs, w s exhibits similarly high values of about 55 m d −1 (Fig. 7f). 60 Inside the OMZs, where remineralization is slower than in oxygenated waters, the detritus residence time is longer and leads to lower w s (Fig. 6l). The subtropical regions, dominated by calcifiers, show a rather homogeneous mean sinking velocity ( w s ≈ 50 m d −1 ) throughout the water column 65 apart from the first few hundred meters and near bottom regions. By contrast, diatom-dominated waters feature a significantly increasing w s with depth, which reaches values of up to approximately 80 m d −1 (Fig. 7f). Diatom-dominated aggregates in surface waters feature a high buoyancy through 70 TEPs and a loose structure which diminishes with continuous remineralization during their descent.
The aggregate properties entering the M 4 AGO scheme are all directly or indirectly measurable. The comparison of simulated and measured aggregate properties is, however, dif-75 ficult, as M 4 AGO depicts mean values of aggregate populations that in situ encompass heterogeneous composition among size spectra. In addition, measurements are often limited to particular aggregate characteristics, while others remain unconstrained within the same data set. We there-80 fore compare M 4 AGO to field and laboratory measurements which examine subsets of the simulated aggregate characteristics.
The modeled aggregate excess densities in diatomdominated regions compare well to former field and labo-85 ratory measurements, where marine aggregates showed excess densities of about 0 to ∼ 10 kg m −3 (Alldredge and Gotschalk, 1988;Ploug et al., 2008;Iversen and Ploug, 2013;Laurenceau-Cornec et al., 2019). In calcifier-dominated regions, aggregate excess densities are about ∼ 35 kg m −3 , 90 which is in the upper range of measured values of 2.1 to 41 kg m −3 (Engel et al., 2009). The increased excess density of CaCO 3 -dominated aggregates of about 100 to 150 kg m −3 in about 1000 m depth compares well with that of observed fecal pellets egested by coccolith-consuming zooplankton 95 (White et al., 2018). These fecal pellets also show similar mean sinking velocities to our modeled w s ≈ 50 m d −1 . The change of the excess density is linked to the increasing d f of aggregates with depth that is in qualitative agreement with present conceptual understanding (Mari et al., 2017). 100 The increasing d f depicts the expected continuous repacking of aggregates and the zooplankton-mediated compaction in fecal pellets. The latter particularly takes place in the upper few hundred meters of the ocean and is a major pathway of coccolith transport (De La Rocha and Passow, 2007;Honjo, 105 1976).
The general difference in typical size between diatomrich aggregates and CaCO 3 shell-enriched aggregates compares well to observations that also showed smaller CaCO 3dominated aggregates (Biermann and Engel, 2010). The de-110 Figure 6. Modeled marine aggregate properties at (a-f) 100 m and (g-l) 960 m depth. Mathematical symbols are as follows: ρ p -mean primary-particle density; d p -mean primary particle diameter; d f -microstructure (fractal dimension); ρ f V -volume-weighted mean excess aggregate density; d max -maximum aggregate diameter; w s -concentration-weighted mean sinking velocity of aggregates. Note that w s comprises the full range of many micrometer-sized to rare, large aggregates with low (O(1 m d −1 )) and high (O(> 100 m d −1 )) sinking velocities. For the mean stickiness, α , volume-weighted mean porosity, φ V , and the number distribution slope, b, see Fig. A1. creasing d max with depth is in qualitative agreement with observed vertically decreasing aggregate mean diameters (De La Rocha and Passow, 2007). Decreasing maximum aggregate sizes with depth and reduced organic matter content also agree qualitatively well with experiments of Hamm (2002) 5 and Passow and De La Rocha (2006). Both studies showed a significant decrease in aggregate size with increasing mineral components. In their experimental setups, it remains elusive if a certain threshold of carrying capacity of POM was reached (Passow and De La Rocha, 2006) or if the bal-10 ance between adhesive forces within the aggregates and the sinking-induced shear forces (Adler, 1979;Brakalov, 1987) was shifted towards smaller aggregates. It is likely that both effects happen at the same time, since natural polymers possess stronger adhesive surface properties than biogenic min-15 erals (Eisma, 1986). In M 4 AGO, only the compaction towards higher fractal dimensions through lower internal binding forces, expressed as stickiness of the primary particles, is represented. Compaction can coincide with an increasing number of binding links in aggregates, which can lower 20 the overall susceptibility of aggregates to shear stress. In M 4 AGO, we disregard this effect and keep Re crit globally constant.
In summary, the resulting mean sinking velocity in M 4 AGO is of same order of magnitude as that found in 25 observations (Alldredge and Gotschalk, 1988;White et al., 2018;Villa-Alfagame et al., 2016). We note, however, that mean sinking velocities estimated from observations potentially overestimate w s , as they (i) are methodologically constrained to particles larger than a lower detection limit, 30 which is typically much larger than primary-particle size, and (ii) depend on often uncertain size-to-mass relationships. We emphasize further that our modeled w s embraces numerous slowly sinking aggregates of primary-particle size (w s ≈ O(1 m d −1 )) as well as rare, but large, fast sinking aggre-35 gates (w s ≈ O(1000 m d −1 )). This range hence encompasses observations for single cells and coccoliths up to large marine snow aggregates and fecal pellets (Miklasz and Denny, 2010;Alldredge and Gotschalk, 1988;Biermann and Engel, 2010). Generally, the spatio-temporal variability in w s in M 4 AGO differs significantly from the simple underlying as-

Global pattern of transfer efficiency
The transfer efficiency of POC from z 0 = 100 m to depth z in the ocean, T eff,POC,z = POC flux at depth z POC flux at depth 100 m , for z > z 0 , provides an estimate on the fraction of exported POC that 15 reaches a particular depth and is determined by sinking velocity and remineralization. In the Martin curve (Eq. 1), the slope constant, β, prescribes the transfer efficiency, to a particular depth and leads to an almost homogeneous 20 global transfer efficiency of T eff,POC,z (Martin) ≈ 0.12 to about 1000 m depth in our standard run (Fig. 8a). The lower remineralization rates in sub-anoxic or even anoxic OMZs com-pared to oxygenated regions lead to higher transfer efficiencies, visible in the equatorial eastern Pacific Ocean, 25 the northern Indian Ocean, CE1 and upwelling regions off the coast of Peru and Africa. Apart from these regions, the standard run features only little variability, as expected from the relationship between the Martin curve slope parameter and the transfer efficiency. The remaining variability is related to 30 ocean currents and spatially variable turbulent mixing. By contrast, the transfer efficiency in M 4 AGO exhibits a distinct global pattern and possesses higher efficiency in highlatitude and upwelling regions compared to the subtropical gyres where low T eff,POC,z appears. Similar to the standard 35 run, the OMZ regions feature high transfer efficiencies in M 4 AGO. Since local Martin curves provide meaningful information on the attenuation of POC fluxes with depth, we analyzed the effective slope parameter, β , for both the standard and the M 4 AGO run. We fitted the Martin curve (Eq. 1), 40 to modeled POC fluxes below z 0 to estimate β . As expected, the standard run shows little spatial variability apart from the OMZs and features a global, area-weighted mean slope of β A,Martin ≈ 0.82. In agreement with the transfer efficiency pattern, M 4 AGO possesses a strong latitudinal pattern of the 45 effective slope β that varies between β ≈ 0.59 and 0.67 in high latitudes (Antarctic zone -AAZ; North Pacific -NP; subantarctic zone -SAZ) and ≈ 0.30 and 0.60 in OMZ regions and with a maximum β ≈ 1.31 in the subtropical Pacific gyres (Fig. B1). allows for a more reliable constrain on POC transfer efficiency (Usbeck et al., 2003;Weber et al., 2016). The comparison of the standard and M 4 AGO run reveals the inherent inability of the Martin approach to capture the latitudinal variability (Fig. 8a-c). M 4 AGO agrees qualitatively and quanti-5 tatively well with the reconstructed transfer efficiency pattern of Weber et al. (2016). The overestimation of the transfer efficiency in the equatorial tropical Pacific (ETP) by both the standard and M 4 AGO run is due to the models overestimation of OMZs extensions (Bopp et al., 2013). The large OMZ 10 causes diminished remineralization and, hence, reduced attenuation of POC fluxes. In general, however, the increased transfer efficiency associated to OMZs is in agreement with observations that suggest lower flux attenuation and, hence, small β (Roullier et al., 2014;Löscher et al., 2016;Le 15 Moigne et al., 2017). Lower remineralization in OMZs keeps the POC-to-ballasting-minerals ratio higher when compared to oxygenated waters. Thus, the lower remineralization rates are potentially accompanied by lower sinking velocities, which provides a positive feedback loop on the OMZ evo-20 lution. The larger the vertical extent of the OMZ becomes, the longer the retention time becomes through higher POM aggregate content and lower sinking velocities. Eventually, OMZ evolution is balanced by the oxygen supply through mixing and transport processes. The global pattern of β in 25 M 4 AGO agrees well with the geographical range and pattern found by Buesseler and Boyd (2009) and suggested by Marsay et al. (2015). However, M 4 AGO shows lower max-imum values of β ≈ 1.31 compared to β ≈ 1.9 suggested by Marsay et al. (2015). Globally averaged, M 4 AGO ex-30 hibits an effective slope parameter of β A,M 4 AGO ≈ 0.75, which is lower than the slope of β = 0.86, originally published by Martin et al. (1987) based on local observations. Noticeably, however, the trend of smaller to larger β values from the nearshore to open waters off the Pacific US coast 35 is in agreement with the elusive trend found by Martin et al. (1987). Overall, M 4 AGO clearly improves the representation of the POC transfer efficiency pattern compared to the standard Martin approach. 40 remineralization to the transfer efficiency pattern

Contributions of w s and temperature-dependent
The attenuation of POC fluxes is primarily regulated by the sinking velocity and total remineralization rate of POC. Depth-dependent, sheared lateral transport and vertical water motion can additionally affect the vertical distribution of 45 particulate matter and thus vertical fluxes. We neglect these processes in the following, since (i) the timescale of sinking from one layer to the next layer below is typically shorter than the lateral transport at grid resolutions used in our model runs, and (ii) sinking velocity is faster than typical vertical 50 motions represented by models with this grid resolution. The remineralization length scale (RLS), z * POC , is given by the local ratio of w s to remineralization (R POC,remin ), We here focus on showcasing the temperature effect on remineralization and calculated the remineralization length scales, z * POC , without oxygen limitation of remineralization, which cancels out for equal O 2 concentrations. Values smaller than 1 imply stronger POC flux attenuation in M 4 AGO than in the standard run. For reference, the RLS in the standard run is given in the top left. The standard RLS increases due to increasing sinking velocity, w * s,POC , with depth (shown bottom left). Contour lines provide the Q 10 factor temperature-dependent remineralization rates with r POC at reference temperature T ref in M 4 AGO. In the standard run, the aerobic rate is globally constant (0.026 d −1 ). (b) Relative contributions of sinking, RC w s , and remineralization, RC remin , to difference between standard and M 4 AGO run. Contour lines provide the relative contribution of remineralization.
and defines the vertical distance in which POC would decay to 1/e (≈ 37 %) of its initial value, the POC e-folding depth (Cram et al., 2018). The RLS can be locally calculated and enables (i) exploring the local effects of remineralization and sinking velocity on the attenuation of POC fluxes and (ii) better understanding the depth-integrated information provided by the transfer efficiency or Martin's effective slope parameter, β . The RLSs in surface waters and the upper mesopelagic 10 zone of subtropical and equatorial regions are shorter in the M 4 AGO run than in the standard run by more than a factor of 2 (Fig. 9a). This higher turnover causes the lower p ratio in M 4 AGO compared to the standard run in these regions (Fig. 4b). In the mesopelagic zone, the RLSs in M 4 AGO are 15 similar or pronounced longer and decrease again in deeper regions compared to the standard run (Fig. 9a). The longer RLSs in the mesopelagic zone of the high-latitude ocean are the reason for the higher transfer efficiency of M 4 AGO compared to the standard run. In order to analyze which of the 20 two processes, sinking or remineralization, is of primary importance for the change in the RLS and thus the transfer efficiency, we define their relative contributions to the difference between M 4 AGO and the standard run as where the processes, P i = { w s , R POC,remin }, refer to Eq. (35), ∂ P i is the partial derivative of z * POC with respect to the process P i (applying the standard run rates), and P i,x is the difference between the value in the M 4 AGO run and the standard run at spatial point x. M 4 AGO possesses gen-30 erally higher remineralization rates than the standard run, which would increase the flux attenuation compared to the standard run. The temperature-dependent remineralization in M 4 AGO shows lower rates in the cold waters of the high latitudes than in the equatorial and subtropical regions (Fig. 9a). 35 Within M 4 AGO, this pattern enhances the RLSs, and thus transfer efficiency, in the high latitudes compared to the equatorial regions, which is in agreement with Marsay et al. (2015) and Cram et al. (2018). Compared to the standard run, higher w s overcompensates the effect of intensified 40 remineralization below the thermocline and leads to a longer RLS (Fig. 9b) in the mesopelagic zone, which is particularly true in diatom-dominated regions. In CaCO 3 -dominated subtropical gyres, w s falls below the sinking velocity of the standard model in regions deeper than 3000 m and thus con-45 tributes further to the stronger flux attenuation compared to the standard run. The higher w s in M 4 AGO turns out to dominate over remineralization and increases the RLS, and hence the transfer efficiency, in the mesopelagic zone of the high latitudes. In summary, the temperature dependence of 50 remineralization in M 4 AGO induces a latitudinal pattern of longer RLSs, and thus higher transfer efficiency, in high latitudes, which is further amplified by high sinking velocities of diatom-dominated aggregates in the mesopelagic zone. The analysis of the RLS changes compared to the standard run emphasizes the role of sinking velocity for the longer RLS in the mesopelagic zone in high latitudes and thus for the enhanced transfer efficiency. We therefore aim to better 60 understand the underlying factors that control the mean sinking velocity. We suggest in the following that the size of pri-mary particles might be as important as the density of the ballasting material for defining the sinking velocity of aggregates.

Impact of mineral size and ballasting effect on
M 4 AGO allows for assessing the contributions of aggregate properties and molecular viscosity to w s . We define the relative contributions of the modeled particle properties and the molecular viscosity that control w s at particular depth z based on a first-order approach: where and R X i ,z = 100, but R X i ,z = 100, and X i,z is the global mean of the contributing property at depth z. The relative contributions provide information about the main driving factors for the local w s when compared to w s of global average aggregates. For example, what is the 15 percentage-wise contribution to the local sinking velocity by local primary-particle density compared to the global mean primary-particle density? By neglecting the higherorder terms, we provide only qualitative insights into the role of the different aggregate properties and molecular viscosity 20 in w s .
CaCO 3 acts as a strong positive, w s -increasing factor for sinking velocities in the equatorial and subtropical regions due to its high primary-particle density (Fig. 10a). As an exogenous factor, the low molecular viscosity in warm regions 25 contributes positively to sinking velocity (Fig. 10d). Additionally, in our model, d f increases particle excess density and thus w s in CaCO 3 -rich areas (Fig. 10c). In the Southern Ocean, North Pacific, North Atlantic and upwelling regions, particularly the large d p in diatom-dominated regions en- 30 hances w s (Fig. 10b). According to our model, d max of aggregates and the number distribution slope, b, contribute positively to w s in the high latitudes, except for the Arctic Ocean, where mineral material strongly affects the aggregate properties. The pattern of highly variable endogenous and ex- 35 ogenous controls on w s emerges particularly in the upper ocean, while in deeper regions of about 1000 m depth, the aggregate properties become more and more homogeneous, except for d p and ρ p and in OMZs (Fig. 10g-l). With homogenization of most aggregate properties with depth, the 40 contrasting relative contributions of dense mean primary particles and large mean primary particles to w s become even more pronounced. Denser, small CaCO 3 particles contribute to w s in the subtropical, oligotrophic regions of the oceans, while comparably less dense but larger opal frustules lead 45 to high w s in nutrient-rich upwelling and high-latitude regions. Even though our formulation for w s is highly nonlinear, we expect the qualitative pattern of the relative contributions to w s , particularly through d p and ρ p , to be coherent. In sum, we therefore emphasize the potential role 50 of primary-particle size, in particular that of diatom frustules, in determining w s and thus POC fluxes.
The effects of microstructure on in situ w s are not well studied, and d f of aggregates is weakly constrained. It therefore warrants further investigation, especially if a spatially 55 variable d f effect occurs on in situ w s in the upper ocean.
Vertically varying d f has been suggested (Mari et al., 2017), and given the general importance of the microstructure for sinking velocity, a deeper understanding of the factors and processes influencing d f is hence highly desirable. 60 Microstructure d f directly prescribes the aggregate number distribution slope, b, in M 4 AGO (see Fig. A1c and f for a  map of b). Observations show a higher variability in b ≈ 2 to 5 (Guidi et al., 2009) than in M 4 AGO, where b varies between ≈ 3.19 and 3.76 and thus likely underestimates the 65 spatial variability and relative contribution of b to w s . At present, M 4 AGO is limited to the steady-state size distribution that represents the characteristic processes of aggregation and fragmentation in the system. Explicit modeling of the dynamics of the aggregate size spectrum would be re-70 quired to cover the variability in measured b, which would further enhance the variability in w s . We briefly discuss this current model limitation in Sect. 3.10.
Previous studies support our finding that CaCO 3 acts as a strong ballasting agent (e.g., Francois et al., 2002;Balch 75 et al., 2010;Cram et al., 2018). By contrast, opal density of the hollow silicate structures has a smaller effect on ρ p (see Fig. 7a). Instead, silicate frustule size of diatoms significantly affects w s (see Figs. 10b and h and 7b). This contrasts the assumption of opal acting solely as ballasting 80 material (Cram et al., 2018) and is congruent with Francois et al. (2002), who noticed that factors other than particle density likely play a role. As indicated by Eqs. (5) and (6), primary-particle size affects the excess density and porosity of aggregates, which have decisive effects on sink-85 ing velocity (Laurenceau-Cornec et al., 2019). Deciphering the size effect of primary particles on in situ w s and fluxes might be particularly challenging in regions with a diverse size structure of the phytoplankton community. Oligotrophic regions typically harbor a narrower size distribution 90 of phytoplankton (see, e.g., size ranges and standing stocks in Kostadinov et al., 2016), which may produce more homogeneous aggregates and thus more predictable POC fluxes, as found by Guidi et al. (2016), who correlated POC fluxes to the oligotrophic phytoplankton community. In addition, 95 interannual variability in the dominant size of primary producers has been suggested to drive the interannual change in export fluxes (Boyd and Newton, 1995). In contrast to oligotrophic phytoplankton communities, diatom-dominated communities feature a higher size diversity (Tréguer et al.,100 2018) and different morphologies, which both affect the sinking velocity of aggregates (Laurenceau-Cornec et al., 2015;Bach et al., 2016). Phytoplankton community size structure and morphology thus introduce higher variability to sinking velocity, which complicates attribution of cell size and mor-105 phology effects on POC fluxes. The higher size and morphology variability in diatom-dominated phytoplankton commu- Figure 10. Qualitative first-order relative contributions of marine aggregate properties to the change of local mean sinking velocity compared to a global mean aggregate spectrum at depth z. (a-f) 100 m; (g-l) 960 m. Mathematical symbols are as follows: ρ p -mean primary-particle density; d p -mean primary-particle diameter; d f -fractal dimension; µ -dynamic molecular viscosity; d max -maximum aggregate diameter; b -aggregate number distribution slope. nities, together with variable cellular silicate-to-carbon ratios (Brzezinski, 1985), likely explains the poor correlations of opal to POC fluxes (e.g., found by Francois et al., 2002). Similarly, the correlation between POC fluxes and, by number and size, more variable foraminiferal CaCO 3 is weaker 5 than for coccolith fluxes (De La Rocha and Passow, 2007). We therefore suggest to factor in, or even focus on, frustule sizes and morphology as explanatory variables for POC fluxes when carrying out mineral ballasting and rain ratio studies in diatom-dominated regions. 10 3.6 Regional fluxes and rain ratios Biogenic minerals, dust particles and detritus are tied together in M 4 AGO and define the composition, microstructure and thus sinking velocity of aggregates. In turn, the tracers are independently remineralized or dissolved. Both pro- 15 cesses affect w s and thus the e-folding depth of tracers. For example, the remineralization of POC increases w s of diatom-dominated aggregates and thus increases the dissolution length scale of opal, which we refer to as the "opal RLS". The combined sinking thus (i) affects the RLS of POC and 20 opal (and CaCO 3 when dissolution takes place) and (ii) couples the timing of mineral and POC fluxes at depth. The Martin curve concept can be applied to represent the effective attenuation of opal fluxes. The effective Martin curve slope for opal, β opal , exhibits similar regional variabil-25 ity in the M 4 AGO run compared to the standard run (Fig. 11). However, the equatorial tropical Pacific (ETP) region shows higher spatial variability. The low POC remineralization in OMZs leads to small w s of aggregates and decreases the opal RLS and thus increases β opal . Similar to POC fluxes, 30 opal fluxes exhibit shorter opal RLSs in the surface waters, while they exceed the standard RLSs in the mesopelagic zone and below (not shown). Particularly in the deep waters of the high latitudes, the opal RLSs are longer than in the standard run. This is due to the aggregates' higher sinking veloc-35 Figure 11. Area-weighted mean effective power function slopes for opal fluxes, β opal A , in the regions defined in Fig. 8. Calculated for opal fluxes below 100 m. Error bars show regional standard deviation.
ity than the 25 m d −1 opal sinking speed in the standard run and partially to the reformulation of temperature-dependent opal dissolution. In the subtropical regions, opal is remineralized faster below 2000 m and fluxes are generally small (see Fig. 5a, b). 5 In sum, M 4 AGO estimates that globally ∼ 1.03 Gt Si yr −1 (∼ 1.04 Gt Si yr −1 in the standard run) reaches the seafloor, which is in good agreement with present estimates of about 1 Gt Si yr −1 (Tréguer, 2002). Generally, the modeled β opal is well within the range of observations with estimates of 10 β opal = 0.22 ± 0.53 (Boyd et al., 2017). In M 4 AGO, dissolution of opal has the implication that ballasting mineral ratios between opal and CaCO 3 shift towards higher importance of CaCO 3 with depth. Since the dissolution of opal is temperature-dependent, the preservation efficiency for opal-15 to-CaCO 3 ratios is thus lowest in regions with strong vertical temperature gradients.
The direct coupling of POC and biogenic mineral fluxes through marine aggregates potentially enhances the models' ability to represent rain ratios in space and time. We there-20 fore aggregated the sediment trap data set of Mouw et al. (2016a, b) with 15 792 individual POC flux measurements at 673 unique locations to a monthly climatology of POC, particulate inorganic carbon (PIC) and silicate fluxes. We accounted for the time spans of sediment trap deployment rang- 25 ing from hours to years by time-weighting with a minimum weight of 1 d per month in cases of few hours of measurement and for each covered month the full monthly weight in cases of a yearly measurement. We compared the climatologies to modeled monthly mean flux ratios at the stations 30 at their respective depths. The M 4 AGO run represents the POC/PIC flux ratio equally as well as the standard run (not shown). For POC/Si fluxes, the scatter around the 1 : 1 line is reduced in M 4 AGO compared to the standard run (Fig. 12).
However, M 4 AGO introduces a bias in regions deeper than 35 ∼ 4000 m towards smaller POC/Si ratios. The RLSs in the deep ocean are hence too short in M 4 AGO. While the scatter of the POC/Si ratio reduces in M 4 AGO, the overall variability in the flux ratios becomes compressed compared to the measured variability. This compression might be caused by 40 several factors. First, we assume a tight connection of aggregate components and ignore heterogeneous composition among local particle-size spectra which would potentially cause different sinking velocities among the different components. Second, HAMOCC ignores changing diatom silici-45 fication caused by temperature (see Ragueneau et al., 2006, and references therein) and by seasonal changes in nutrient availability (Assmy et al., 2013). Generally, deficits exist in both the model runs and sediment trap data. The lack of resolving small-scale variability through eddies and their role 50 in shaping the phytoplankton community, as well as general timing, internal variability, and spatial current shifts, limit the global models' ability to represent local features. Measurement limitations of sediment traps have the potential to additionally increase the mismatch (see, e.g., Usbeck et al., 55 2003, and references therein).

Evaluation of biogeochemical tracer distributions
We evaluate the simulated spatial distribution of nutrients, oxygen and alkalinity by comparison to gridded observation climatologies. Silicate, phosphate, nitrate and oxygen 60 are compared to the World Ocean Atlas (WOA; Boyer et al., 2013;Garcia et al., 2014b). Alkalinity, A T , is compared to the Global Ocean Data Analysis Project (GLODAPv2) climatology (Lauvset et al., 2016). We present the results in Taylor diagrams (Taylor, 2001), which aggregate correlation, 65 root-mean-square deviation (RMSD), and standard deviation of simulated and observed variables into distances between model results and the reference observations ( Fig. 13a-d).
Here, we use spatial grid-cell-wise statistics derived from temporally averaged (100 year mean) fields for which we 70 interpolated the observational data to the model grid. To be able to show all parameters in one plot, we derive normalized standard deviations and RMSDs. We evaluate the Atlantic and the Pacific Ocean separately to account for their different hydrographical features. For example, the Atlantic Ocean is 75 characterized by ventilation through deep water formation in the high latitudes, a feature that does not exist in the Pacific. Furthermore, we select four depth levels: surface waters (6 m depth); the two depths that determine the transfer efficiency, 100 m and approximately 1000 m; and an intermediate depth 80 in the mesopelagic zone, 362 m. Generally, the two model runs represent nutrients, silicate and oxygen distributions equally well at these depth ( Fig. 13a-d); i.e., differences to observations are larger than among the two models. This is expected, as tracer distribu-85 tions in the model are predominantly determined by the flow field and both simulations use identical physical model conditions. The Taylor diagram captures matches of spatial pattern such as location of fronts, extension of gyres and the location of water masses in an aggregated form. The ocean 90  Fig. 5a, b). Correlation coefficient, r, and significance value, p, are given if r > 0.4. model represents such features realistically (Jungclaus et al., 2013), but we cannot expect a perfect match of the circulation pattern with in situ conditions for multiple reasons. The climatological, simplistic atmospheric forcing damps in particular the interannual variability in ocean currents and can con-5 tribute to spatial shifts of, for example, frontal regions with respect to the observed climatology. In addition, the coarse horizontal resolution of the model and atmospheric forcing affects features such as upwelling strength (Milinski et al., 2016) and location. In turn, the spatio-temporal interpolation 10 of data, necessary through scarcity of observations, limits the resolution of nutrient variability and introduces uncertainty, particularly in deeper regions, where the density of observations is lower than in surface waters. As a consequence, a mismatch between modeled biogeochemical tracers and the 15 climatological mean of observations is expected.
Tracers undergoing less complex biogeochemical cycling, such as phosphate or silicate, reflect the quality of the flow field more directly than tracers such as nitrate, oxygen or alkalinity (see also England and Maier-Reimer, 2001). The 20 latter tracers are therefore generally more prone to biogeochemical model reformulation, and larger differences can be expected between M 4 AGO and the standard run. In addition, M 4 AGO links the cycles of phosphate and silicate closer to nitrate and oxygen through common sinking of particulate 25 matter. Phosphate and silicate therefore experience an additional biogeochemical cycle imprint in M 4 AGO compared to the standard run. In surface waters of the Pacific, this led, for example, to a slight improvement of the phosphate correlation with WOA observations, while it increased the 30 normalized standard deviation compared to the standard run (Fig. 13a). The improvement is related to the better representation of surface phosphate concentrations off the coast of Peru and the western coast of North America (see Fig. 13eg). 35 In both simulations, biogeochemical tracer distributions in the Pacific are strongly influenced by the OMZ in the eastern boundary upwelling region. An exception is oxygen in surface waters, which is primarily determined by gas exchange processes. As a common feature of all state-of-the-art 40 Figure 13. (a-d) Taylor diagram for tracers (oxygen -O 2 ; nitrate -NO − 3 ; phosphate -PO 3− 4 ; silicate -Si) in comparison to the World Ocean Atlas data (Boyer et al., 2013;Garcia et al., 2014b) and, in the case of total alkalinity, A T , to GLODAPv2 data (Lauvset et al., 2016). global ocean biogeochemistry models, OMZs are too large and result from nutrient trapping through a sluggish circulation and associated insufficient supply of oxygen (Dietze and Loeptien, 2013). The sluggish circulation and mixing are likely associated with underrepresented equatorial cur-5 rents, particularly the equatorial intermediate current system and equatorial deep jets (Brandt et al., 2012;Shigemitsu et al., 2017). The OMZ shape in the eastern Pacific Ocean, though, changes between the runs. In M 4 AGO, the water column above 800 m is more oxygenated and below 800 m 10 is less oxygenated than in the standard run. The change in OMZ shape imprints on A T , as anaerobic and aerobic remineralization processes change alkalinity in a different man-ner; i.e., denitrification and sulfate reduction increase alkalinity, whereas aerobic remineralization decreases alkalinity. 15 While A T and silicate are generally worse than other tracers when compared to observations, M 4 AGO improves the representation of A T in the Pacific at 360 m compared to the standard run. This is visible through the higher correlation (∼ 0.45 compared to 0.4) and lower RMSD of ∼ 2.3 com-20 pared to ∼ 3.75 in the standard run ( Fig. 13c; 362 m depth). The increased RMSD of silicate in surface waters and 100 m in M 4 AGO are associated to the eastern equatorial and upwelling regions in the Pacific. The lower remineralization in OMZs not only decreases the sinking velocity of detritus but 25 also increases the retention time of the tightly coupled opal in M 4 AGO and thus enhances the opal dissolution and silicate concentration. In general, the tighter coupling of silicate to nutrients through common sinking in aggregates leads to a higher sensitivity of silicate to ocean circulation and temperature deficiencies in M 4 AGO than in the standard run. The high RMSD and low correlation for silicate compared to observations, particularly in the euphotic zone, can have additional implications for M 4 AGO, since the silicate distribution in the euphotic zone directly affects opal production and thus the sinking velocity of aggregates, transfer efficiency and nu-10 trient distributions. The tighter coupling of silicate to other biogeochemical cycles in M 4 AGO therefore warrants further future investigation.

Regional CO 2 fluxes
In a 100-year climatological steady state, both model runs 15 show the Southern Hemisphere ocean acting as a net source of CO 2 to the atmosphere, while the Northern Hemisphere acts as a net sink (Fig. 14a). Consequently, net oceanic CO 2 transport across the Equator from the Northern to the Southern Hemisphere exists. In the simulation with M 4 AGO, a 20 stronger CO 2 uptake compared to the standard run occurs in the region between 60 and 45 • S, which coincides with deeper transfer of POM reflected by an increased transfer efficiency (Fig. 8). In both model runs, the tropical regions outgas CO 2 to the atmosphere. As is visible from the decreas- 25 ing difference between the cumulative zonal CO 2 fluxes, M 4 AGO exhibits stronger outgassing in the tropical region, particularly in the Northern Hemisphere (Fig. 14a). Even though a small residual of zonally integrated fluxes of about 0.02 Gt C yr −1 between the two model runs remains, the small residual CO 2 flux imbalance of 0.07 and 0.05 Gt C yr −1 in M 4 AGO and the standard run, respectively, indicate wellspun-up model runs. In general, the latitudinal zonal CO 2 fluxes and the cross-equatorial southward oceanic CO 2 transport agree qualitatively and quantitatively well with former 35 forward-integrated models for preindustrial conditions (e.g., Sarmiento et al., 2000;Gloor et al., 2003;Mikaloff Fletcher et al., 2007).
To study the effect of the changed transfer efficiency in M 4 AGO on the CO 2 uptake without indirect climate effects, 40 we perform model runs in which we linearly increase the atmospheric CO 2 concentration with a yearly increment from ∼ 285 to 400 ppm within 150 years. The transient forcing does not feed back onto the physical climate state. We therefore emphasize that no climate feedback on, for example, 45 ocean temperature, stratification, and thus primary production and remineralization occurs in these model runs.
During the atmospheric CO 2 increase, both model runs show similar global annual oceanic CO 2 uptake (Fig. 14b). About 181-185 Gt C is taken up by the global ocean through-50 out the 150-year period. The discrepancy of ≈ 4 Gt C between the two runs arises to a large extent from the resid- Figure 14. (a) Climatological cumulative zonal CO 2 flux in the standard and the M 4 AGO run (from south to north). Generally, negative fluxes represent net CO 2 uptake by the ocean. Note the second y axis for the difference between the two runs (also in ck). (b) Global annual net CO 2 flux under linearly increasing atmospheric CO 2 concentration. (c-k) Regional cumulative CO 2 fluxes under increasing atmospheric CO 2 (regions as defined in Fig. 8).
The general trends in the regional cumulative fluxes re-55 main the same between the two model runs (Fig. 14c-k), and most regions act as net sinks of CO 2 . Only the upwelling regions of the equatorial tropical Atlantic (ETA) and Pacific (ETP) are net sources of CO 2 to the atmosphere at the end of the 150-year period. The magnitude of the regional cumula-60 tive fluxes, however, varies among the two model runs.
Since we can rule out temperature-driven effects on the differences in CO 2 uptake when comparing the two climatologically forced runs, a different state and/or changes in A T and dissolved inorganic carbon (DIC) concentration deter-65 mine the partial pressure of CO 2 , pCO 2 and thus CO 2 fluxes. In our runs, major changes in A T and DIC can locally occur due to a change in, first, the CaCO 3 -to-detritus rain ratio as a result of varying CaCO 3 and POM production and remineral-J. Maerz et al.: Microcomposition of marine aggregates as co-determinant for global POC fluxes ization, or, second, an increase in the POC transfer efficiency, for example, which reflects a transfer of carbon to deeper waters where CO 2 is withdrawn from immediate exchange with the atmosphere. In addition, such different signals can be transported by the flow field, and remote effects can appear. Quantitatively, differences of regional cumulative CO 2 fluxes larger than 5 Gt C appear in the Antarctic zone (AAZ), the subtropical Pacific (STP), the equatorial tropical Pacific (ETP) and the Indic Ocean (IO; Fig. 14c-k). In the upwelling and high-latitude regions, the pCO 2 is lower in M 4 AGO than in the standard run, which translates to higher oceanic uptake (in AAZ and NP) or less outgassing (in ETP). Qualitatively, this coincides well with the primary production, respective export and the higher transfer efficiencies in these regions (cf. Figs. 4 and 8). In regions of higher transfer efficiency but 15 similar CO 2 fluxes when compared to the standard run, either POC export fluxes are small (e.g., in the Arctic Ocean) or physical processes such as mixing or upwelling dominate over biologically induced CO 2 fluxes (e.g., in the SAZ). In the STP region, where high CaCO 3 -to-detritus rain ratios oc-20 cur (refer to Fig. 5), A T is about 0.2 % lower in M 4 AGO than in the standard run, which enhances pCO 2 by about 2 % and thus decreases the CO 2 uptake.
In summary, under linearly increasing CO 2 in a noninteractive mode, the global CO 2 fluxes remain similar for 25 the standard and the M 4 AGO run. However, regional CO 2 fluxes change, which is potentially linked to the changing pattern of transfer efficiency and the CaCO 3 -to-detritus rain ratio. A detailed study on the feedbacks between the transfer efficiency, represented by M 4 AGO, and a transient climate 30 on oceanic CO 2 uptake is out of the scope of this paper and will be part of future investigations.

Sensitivity analysis
With M 4 AGO, the processes of temperature-dependent remineralization and a number of new model parameters to rep-35 resent variable sinking velocity of marine aggregates were introduced to HAMOCC. We carried out three sensitivity experiments to provide insights into the response of HAMOCC to the uncertainty of selected parameters (see below for the criteria). As target variables, we chose (i) the transfer ef- 40 ficiency, (ii) phosphate as an essential nutrient for primary production, and (iii) silicate as a nutrient and circulationreflecting agent.
Previously in Sect. 3.5, we suggested the size of primary particles, particularly that of opal frustules, to be an important factor for regulating w s and thus the transfer efficiency. We thus study the effect of primary-particle size in the sensitivity experiment S(d p,frustule ) exemplarily for diatom frustules because (i) diatom size is highly variable, (ii) HAMOCC does not explicitly represent algae size 50 classes and (iii) algae body size is expected to decrease with rising water temperatures (Daufresne et al., 2009). We decrease their size by 50 %, from d p,frustule =20 µm to d p,frustule =10 µm. In the second sensitivity experiment, S(d f ), the varying fractal dimension, ranging between 1.6 55 and 2.4 in the M 4 AGO run, is set to a constant, d f = 2, to eliminate its variable effect on b and w s . In the third experiment, S(K P,diaz ), we showcase the sensitivity of phosphate concentration in the subtropical gyres to diazotrophs' standing stock, which we noticed during the process of tun-60 ing the M 4 AGO run. Here, we increase the half-saturation constant for phosphate uptake by diazotrophs from 0.05 to 0.10 µmol L −1 . All sensitivity runs are in a quasi-steady state, as reflected by surface properties that are in steady state.
Decreasing the opal frustule size by 50 % compared to 65 the M 4 AGO run reduces w s and thus the RLS in diatomdominated regions. As a consequence, the transfer efficiency in silicifier-dominated regions is lower in S(d p,frustule ) than in M 4 AGO (Fig. 15g). Accordingly, opal dissolves closer to surface waters and the silicate concentration increases com-70 pared to the M 4 AGO run in silicifier-dominated regions, particularly in and downstream of coastal upwelling regions (Fig. 15d). In the subtropical gyres, where silicate is diminished, it further decreases in S(d p,frustule ). The higher remineralization in the euphotic zone in upwelling regions in-75 creases the phosphate concentrations downstream in the subtropical gyres (Fig. 15a). In sum, the sensitivity to d p,frustule underpins the potential importance of primary-particle size for w s and thus for biogeochemical cycling.
Generally, eco-physiological responses of primary pro-80 ducers are typically neglected in ESM-type models such as HAMOCC. Under the premise that the size of primary particles is of importance to represent w s as depicted by M 4 AGO, changes of the phytoplankton size structure with ongoing ocean warming could affect RLSs, transfer effi-85 ciency and thus the biological carbon pump. Indeed, body size decreases with increasing water temperature (Daufresne et al., 2009). The increasing water temperature has been suggested to shift eco-physiological regions polewards by rates of about 22 to 36 km per decade (Lefort et al., 2015, and ref-90 erences therein). S(d p,frustule ) shows that decreasing primaryparticle size reduces the transfer efficiency and thus likely weakens the biological carbon pump. Yet, the net effect of such eco-physiological responses and adaptability of primary producers on transfer efficiency and CO 2 uptake under ongo-95 ing climate change remains elusive and thus demands further future investigation.
In the sensitivity study S(d f ), d f in surface waters is increased in diatom-dominated waters and reduced in calcifierdominated waters compared to the M 4 AGO run (cf. Fig. 6c). 100 In the M 4 AGO run, aggregates experience rapid compaction within the depth of about 150 to 250 m in diatom-dominated regions, d f increases rapidly (see Fig. 7c) and, hence, so does w s . By contrast, w s is only enhanced in the productive surface waters of the SAZ, upwelling regions and 105 the associated OMZs due to the larger d f and the dependent smaller b in S(d f ). The lower w s in large parts of the mesopelagic zone and below leads to a generally de- Figure 15. Sensitivity of (a-c) surface phosphate, (d-f) silicate concentrations, and (g) POC transfer efficiency to diatom frustule size, d p,frustule (S(d p,frustule ), modified from 20 to 10 µm), the half-saturation constant for phosphate uptake by diazotrophs (S(K P,diaz ), modified from 0.05 to 0.1 µmol L −1 ), and fractal dimension, d f (S(d f ), from variable to constant d f = 2). Error bars in (g) represent regional standard deviation and, for Weber et al. (2016), the uncertainty for the reconstruction of the regional transfer efficiency. creased transfer efficiency in S(d f ) (Fig. 15g). Consequently, the RLSs are shorter and enhance the silicate concentrations in surface waters in diatom-dominated regions (Fig. 15e). In the subtropical gyres, the silicate diminishes even further, which is likely related to the enhanced bulk phytoplankton 5 growth through higher phosphate concentrations. Phosphate is remineralized more rapidly than silicate and thus can be transported downstream the equatorial current into the subtropical gyres, where it is enhanced in S(d f ) compared to M 4 AGO (Fig. 15b). As a consequence of the potential sen-10 sitivity of biogeochemical cycles on d f , more knowledge of processes affecting the aggregate microstructure is required, among them, for example, are compaction through repacking and zooplankton egestion of fecal pellets.
Diazotrophs in HAMOCC can grow independently of ni-15 trate on phosphate. They are regionally confined to warm tropical and subtropical regions (Paulsen et al., 2017(Paulsen et al., , 2018. Diazotrophs modulate the phosphate concentration in the subtropical gyres. The diazotrophs' global primary production increases from about 2.2 Gt C yr −1 in the standard run 20 to ≈ 2.9 Gt C yr −1 in the M 4 AGO run. A slower growth response to phosphate concentrations, enforced by raising the half-saturation constant from 0.05 to 0.10 µmol L −1 in S(K P,diaz ), increases the phosphate concentrations in the subtropical gyres to more than 150 % (Fig. 15c). The primary 25 production through diazotrophs reduces significantly from a former global value of ≈ 2.9 to ≈ 1.9 Gt C yr −1 in S(K P,diaz ) and regionally particularly in the equatorial Panama Basin. By contrast, diazotrophs feature primary production of more than 5 and 7 Gt C yr −1 in S(d p,frustule ) and S(d f ), respectively, 30 due to the increased phosphate concentrations in the subtropical gyres. Diazotrophs in the model only produce organic matter. In S(K P,diaz ), the lower diazotroph primary production in the equatorial upwelling regions (ETA and ETP) thus leads to a shift of the aggregate composition towards higher 35 opal-to-detritus ratios, which slightly increases the transfer efficiency and reduces the available silicate (Fig. 15f). Phosphate previously utilized by diazotrophs in the Panama Basin now partially populates the downstream equatorial current and reaches the subtropical gyres. In comparison 40 to the standard run, this phenomenon is generally intensified in M 4 AGO through the shallower remineralization and lower transfer efficiency in the subtropical gyres (Fig. 8c).
In conclusion, the representation and effects of diazotrophs in HAMOCC, particularly in conjunction with M 4 AGO, re-45 quire further future evaluation.

Current limitations of M 4 AGO
Developing M 4 AGO, we followed a process-oriented approach and explicitly incorporated ballasting and microstructure of aggregates of heterogeneous composition to calculate the mean sinking velocity. Introducing such complexity in 5 ESMs typically comes at the cost of high computational efforts. This is a non-negligible factor for model development which we reduce to a minimum with M 4 AGO. Acknowledging the trade-off between increasing model complexity and computational limitations lets us deploy a number of simplifying assumptions that we critically review in the following.
In the euphotic zone of the oceans, changing phytoplankton community structure, phytoplankton growth, decay and grazing through zooplankton lead to a dynamic supply of various small particles that potentially aggregate and eventually 15 sink and become remineralized. Assuming homogeneous composition of aggregates throughout a dynamic steady-state size distribution is therefore only the first step towards a model representation of marine snow, where local diversity of particles and the processes of aggregation and fragmenta-20 tion are explicitly resolved. By assuming a dynamic steady state for the size distribution, we only represent the characteristic processes within the system and underestimate the variability in the number distribution slope, which ranges between b ≈ 3.19 and 3.76, compared to measurements, where 25 b ≈ 2 to 5 (Guidi et al., 2009). As a consequence, M 4 AGO probably underestimates the variability in sinking velocity, particularly in the upper ocean, where the size distribution dynamically evolves.
In M 4 AGO, TEPs are simplistically considered. In agree-30 ment with Mari et al. (2017), TEPs are microstructureloosening and buoyancy-adding for diatom-dominated aggregates. For CaCO 3 , we only considered the coccolith size which decreases sinking velocity compared to an aggregate composed of intact coccospheres. For simplicity, we 35 assumed that aggregated TEP and intact cells sink as fast as coccoliths and detritus. To decipher the role of TEPs for aggregate formation (Passow, 2002;Mari et al., 2017), our model would require (i) representing TEPs and disintegrating coccospheres and (ii) considering aggregation and fragmen-40 tation processes explicitly to depict the dynamic evolution of marine aggregate distributions. Such a dynamic approach is necessary for, for example, studying short-term events of high POC export hypothesized to be driven by silicate depletion and TEP release by diatoms and their subsequent ag-45 gregation (Martin et al., 2011). More detailed models, such as, for example, the 1-D Lagrangian approach of Jokulsdottir and Archer (2016), can likely provide more insights into aggregate dynamics and can help to further improve the aggregate representation in ESM frameworks. Furthermore, an 50 explicit representation of aggregation and fragmentation enables transient size distributions and would allow for a direct comparison to a growing number of particle distribution slope measurements (e.g., Guidi et al., 2009). 55 We implemented the novel Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M 4 AGO) scheme in HAMOCC to improve the representation of the biological carbon pump in an Earth system model framework. M 4 AGO accounts for the heterogeneity and mi-60 crostructure of aggregates and thus clearly defines measurable statistic aggregate properties in HAMOCC. M 4 AGO links the nutrient and silicate cycle closer together by incorporating opal and other ballasting minerals in aggregate formation and sinking. This lets us introduce a consistent Q 10 65 temperature-dependent dissolution and aerobic remineralization of opal and POC, respectively. In contrast to the standard HAMOCC version, M 4 AGO provides a mechanistic understanding for the recently published global transfer efficiency pattern and represents it 70 well. We identify primary-particle size, particularly of diatom frustules, and the compaction of aggregates with depth as strong driving factors for sinking velocity that, in combination with temperature-dependent remineralization, codetermine the high POC transfer efficiency in high latitudes 75 and upwelling regions. Our model results support previous findings that CaCO 3 with its high density acts as ballasting mineral in calcifier-dominated regions of the ocean. The changed transfer efficiency pattern in combination with the CaCO 3 -to-detritus rain ratio alters regional CO 2 fluxes, while 80 the global uptake remains the same as in the standard run, when atmospheric CO 2 is linearly increased without climate feedbacks.

Summary and conclusions
The highest uncertainties in parametrizing M 4 AGO are with respect to the weakly constrained primary-particle sur-85 face properties and their likely effect on the related microstructure of aggregates. Since sinking velocity and transfer efficiency are highly sensitive to the microstructure of aggregates, gaining insights into the controlling factors and processes for microstructure is desirable. Future model de-90 velopment for the representation of marine aggregates would highly benefit from sub-aggregate-scale measurements of microstructure, adhesive surface properties and primaryparticle composition.
Our findings and the underlying model concept suggest 95 a number of implications. First, the finding that the size of aggregate constituents, particularly of diatom frustules, acts as a potential factor for high sinking velocities, suggests widening the perspective of mineral ballast studies towards a size-and-ballast hypothesis. Further, such an ex-100 tended size-and-ballast hypothesis requires factoring in the different temperature-dependent remineralization and dissolution rates that aggregates experience during their descent. Accounting for cell size and morphology will aid in better assessing the role of the phytoplankton community size struc-105 tures in POC fluxes, particularly in nutrient-rich upwelling regions, where a wide, variable size spectrum of diatoms prevails. Second, the indirect temperature effect on phytoplank-ton cell size (e.g., Daufresne et al., 2009) poses the challenging task of resolving the temperature-adapting cell size structure of the phytoplankton community in global carbon cycle models to depict the potential effect on sinking velocity and thus the biological feedback on rising CO 2 in the atmosphere. 5 In conclusion, M 4 AGO provides well-defined aggregate properties in an ESM framework and can thus serve as a test bed for upscaling aggregate-associated processes to potential global impacts on biogeochemical cycles and, in particular, on the biological carbon pump.
J. Maerz et al.: Microcomposition of marine aggregates as co-determinant for global POC fluxes Appendix A: Additional aggregate properties In addition to the aggregate properties presented in Fig. 6, M 4 AGO involves the aggregate mean stickiness and the number distribution slope that are tightly connected to each other via d f (Fig. A1). In addition, the porosity of aggre-5 gates (Eq. 6) can be deduced from aggregate size, primaryparticle size and d f . Porosity is frequently calculated for aggregates (Alldredge and Gotschalk, 1988;Ploug et al., 2008). Its potential dependency on the microstructure and primary-particle size is, however, seldom covered. We there-10 fore calculated the mean volume-weighted porosity of aggregates, φ V , to provide perspective on how aggregate porosity varies with the aggregate properties shown in Fig. 6 in the global ocean (Fig. A1).
As defined by the attributed primary-particle stickiness, 15 aggregates in diatom-dominated regions show the highest stickiness values in surface waters, since the virtual TEP particles linked to detritus increase the α in M 4 AGO ( Fig. A1). At the bottom of the mesopelagic zone, α is almost homogeneous apart from the OMZ regions, where the lower 20 anaerobic remineralization retains higher detritus concentrations that lead to higher α in conjunction with diatom frustules. The number distribution shows an inverted picture compared to α with the smallest decay slope, b, in diatom-dominated and OMZ regions and the strongest de- 25 cline in calcifier-dominated surface waters (Fig. A1c, f). In M 4 AGO, φ V is tightly connected to α in the euphotic zone (Fig. A1b). With increasing depth, aggregates get compacted, d f increases and the maximum size of aggregates decreases. Both lead to lower φ V , particularly in the diatom-30 dominated regions. In OMZs, however, φ V is high. Generally, the exceptional behavior of aggregate properties in OMZs due to implicitly modeled TEPs is likely overestimated. TEPs possess higher remineralization rate than detritus (Mari et al., 2017), which likely reduces the TEPs oc- 35 currence in deep OMZs and thus their influence on aggregate properties. The effective Martin curve slope β provides a meaningful measure on the attenuation of POC fluxes. This made β a widely used measure to evaluate POC concentration and flux observations (e.g., Lutz et al., 2007;Lam et al., 2011;Marsay 5 et al., 2015). In M 4 AGO, we found a β pattern similar to the inverse of the transfer efficiency ( Fig. B1 compared to Fig. 8, respectively). The smallest β values are found in the models OMZ regions and the shallow Arctic shelf regions.
In the North Pacific and North Atlantic, β features values 10 of about 0.47 to 0.60 and reaches maximum values of about 1.31 in the subtropical gyres. Qualitatively, the pattern thus follows and underpins the previously suggested POC flux attenuation pattern (Marsay et al., 2015;Weber et al., 2016;DeVries and Weber, 2017). By contrast, in our standard run, 15 we found, apart from OMZs, a rather homogeneous β with a global value of β A,Martin ≈ 0.82 that is smaller than the prescribed value of β = 1.00 in oxygen-saturated waters. This discrepancy can be attributed to a number of processes. First, the global β A,Martin is reduced by the lower anaer-20 obic remineralization in OMZs. Second, turbulent diffusion and vertical transport processes represented by MPIOM in addition to sinking contribute to vertical POC concentration profiles and fluxes (Boyd et al., 2019). Third, the artificial numerical diffusion inherent in HAMOCC's implicit upstream 25 scheme for particle sinking contributes to higher mass transport to depth than prescribed by β = 1. This inherent numerical diffusion is, however, implicitly accounted for during the process of tuning β in the model.  The inversely identified transfer efficiency pattern by Weber et al. (2016) provides an estimate on time-integrated climatological POC fluxes. In regions of high seasonal variability in primary production, POC fluxes can undergo strong seasonal 5 variation, and so does transfer efficiency (Lutz et al., 2007). For example, if we assume an average sinking speed of about 25 m d −1 , a pulsed flux at 100 m reaches the depth of 1000 m about a month later and can strongly alter flux and concentration profiles (Lam et al., 2011;Giering et al., 2017). Accord-10 ingly, single measured POC concentration profiles can even show higher concentration at depth than in surface waters, which results in negative β values and thus higher transfer efficiency than 1. In turn, transfer efficiency can be extremely small at the beginning of a phytoplankton bloom. Seasonal 15 transfer efficiency can thus heavily deviate from the climatological state. This seasonal behavior is reflected in M 4 AGO, which shows high transfer efficiency in late autumn and early times of low primary production after the bloom in high latitudes (Fig. C1). Even though the standard run exhibits a sim-20 ilar qualitative pattern, the seasonal amplitude of the transfer efficiency in high latitudes is lower in M 4 AGO compared to the standard run.  nificant discussions on aggregate microphysics. All authors of the paper critically discussed the presented results and contributed by providing valuable feedback during the paper compilation.
Competing interests. The authors declare that they have no conflict of interest. 20 Acknowledgements. The authors thank Thomas Weber for sharing the results on the transfer efficiency and Adrian Burd for the discussion on fractal dimension of aggregates. The authors thank Bo Liu for the internal review and comments on the paper. The authors thank the two anonymous reviewers for their constructive and valu- Financial support. This research has been supported by the Euro- 35 pean Commission H2020 Research Infrastructures (CRESCENDO (grant no. 641816)).
The article processing charges for this open-access publication were covered by the Max Planck Society. 40 Review statement. This paper was edited by Carolin Löscher and reviewed by two anonymous referees.