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 implemented a novel Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M⁴AGO) sinking scheme. M⁴AGO explicitly represents the size, microstructure, heterogeneous composition, density and porosity of aggregates and ties ballasting mineral and particulate organic carbon (POC) fluxes together. Additionally, we incorporated temperature-dependent remineralization of POC. We compare M⁴AGO with the standard HAMOCC version, where POC fluxes follow a Martin curve approach with (i) linearly increasing sinking velocity with depth and (ii) temperature-independent remineralization. Minerals descend separately with a constant speed. In contrast to the standard HAMOCC, M⁴AGO reproduces the latitudinal pattern of POC transfer efficiency, as recently constrained by Weber et al. (2016). High latitudes show transfer efficiencies of $\approx \mathrm{0.25}\pm \mathrm{0.04}$, and the subtropical gyres show lower values of about 0.10±0.03. In addition to temperature as a driving factor for remineralization, diatom frustule size co-determines POC fluxes in silicifier-dominated ocean regions, while calcium carbonate enhances the aggregate excess density and thus sinking velocity in subtropical gyres. Prescribing rising carbon dioxide (CO₂) concentrations in stand-alone runs (without climate feedback), M⁴AGO alters the regional ocean atmosphere CO₂ fluxes compared to the standard model. M⁴AGO exhibits higher CO₂ uptake in the Southern Ocean compared to the standard run, while in subtropical gyres, less CO₂ is taken up. Overall, the global oceanic CO₂ uptake remains the same. With the explicit representation of measurable aggregate properties, M⁴AGO can serve as a test bed for evaluating the impact of aggregate-associated processes on global biogeochemical cycles and, in particular, on the biological carbon pump.

1 Introduction

Marine aggregates transfer biologically bound carbon and nutrients from the sunlit surface waters, the euphotic zone, to the interior of the oceans. While uncertainty with respect to primary production estimates exists, about 4.0 to 11.2 Gt C yr⁻¹ 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 net withdrawal of carbon dioxide (CO₂) 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 biological 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 (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 quantifying 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.

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. 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 carbonate (CaCO₃) 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 capability for ballasting minerals before aggregates disintegrate (Passow, 2004; 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 out of the euphotic zone that reaches a particular depth, e.g., 1000 m (Francois et al., 2002). As CaCO₃ is significantly denser than opal, CaCO₃ is suggested to be a more effective ballasting material for marine aggregates (Balch et al., 2010), implying higher POC transfer efficiency in CaCO₃-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₃-to-opal ratios are found in oligotrophic regions of the mid-latitude subtropical gyres, while 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; 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 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 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 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 understanding 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 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₃ and opal, temperature effects on microbial aerobic and anaerobic remineralization, 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 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 typically formulated to reproduce the Martin curve (Martin et al., 1987) or heuristically describe ballasting of POC with opal and CaCO₃ (e.g., Gehlen et al., 2006; Heinemann et al., 2019), which limits the process-based adaptation of sinking velocities under changing environmental conditions associated with climate change.

As a first step, we develop the Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M⁴AGO) sinking scheme that explicitly represents composition, microstructure, and related properties such as porosity and density of aggregates. We aim to consistently define marine aggregates with their in situ measurable properties in an ESM framework. We implement M⁴AGO in the global HAMburg Ocean Carbon Cycle (HAMOCC) model, which is part of the Max Planck Institute – Earth system model (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 the underlying driving factors for this pattern and (iii) give insights into the impact of M⁴AGO on the global CO₂ 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 concentration 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., 2012; Marsay et al., 2015).

With M⁴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 rates.

2 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; Ilyina et al., 2013; 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 experience 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, 2018). Zooplankton feeds on bulk phytoplankton and releases POM, which enters the common detritus pool. During detritus formation through bulk phytoplankton or zooplankton, opal or CaCO₃ is produced depending on silicate availability. This treatment adequately depicts the spatial distribution of silicifying 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 and evaluated in previous studies; for details, see, e.g., Six and Maier-Reimer (1996), Ilyina et al. (2013), 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 M⁴AGO sinking scheme. A table with the used mathematical symbols can be found in Appendix D, Table .

2.1 HAMOCC standard representation of sinking fluxes and remineralization

The standard version of HAMOCC (Mauritsen et al., 2019) represents sinking fluxes of POC, F_POC, at depth z>z₀, according to the concept of the Martin curve (Martin et al., 1987; Kriest and Oschlies, 2008):

$$\begin{array}{}\text{(1)}& F_{\mathrm{POC}}\left(z\right)=F_{\mathrm{0}}{\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-\mathit{\beta }},\end{array}$$

where F₀ is the POC flux out of the euphotic zone at export depth z₀. For simplicity, the export depth is globally defined as being z₀=100 m in HAMOCC. Above z₀, a constant sinking speed of 3.5 m d⁻¹ is assumed. Below z₀, we assume a linearly increasing mass concentration-weighted mean sinking velocity with depth. The ratio between the remineralization rate of POC, R_POC,remin, and the vertical gradient of the sinking velocity, ${\partial }_z\stackrel{\mathrm{‾}}{w_{\mathrm{s}}}$, determines the POC flux slope, $\mathit{\beta }=R_{\mathrm{POC},\mathrm{remin}}/{\partial }_z\stackrel{\mathrm{‾}}{w_{\mathrm{s}}}$ (Kriest and Oschlies, 2008). In the standard version of HAMOCC, remineralization of POC is temperature-independent and comprehends oxygen-concentration-dependent aerobic remineralization as well as sulfate reduction and denitrification under sub-anoxic and anoxic conditions.

The sinking tracers, opal and CaCO₃, are treated separately from POC and sink with their own, constant sinking 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 conditions below the dynamically emerging lysocline. In the following, we refer to this version of HAMOCC as “standard”.

2.2 The novel M⁴AGO sinking scheme in HAMOCC

Natural waters exhibit a size spectrum of aggregates whose diameter, d, composition and microstructure determine their terminal sinking velocity. In the M⁴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 (following e.g., Kriest and Evans, 1999; Gehlen et al., 2006; Schwinger et al., 2016):

$$\begin{array}{}\text{(2)}& n\left(d\right)=a{d}^{-b}.\end{array}$$

This way, we avoid the computational costs of size-class-based model approaches (Jackson, 1990; Stemmann et al., 2004; Sherwood et al., 2018).

The local concentration-weighted mean sinking velocity, 〈w_s〉, in M⁴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 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〉

$$\begin{array}{}\text{(3)}& \langle w_{\mathrm{s}}\rangle ={\displaystyle \frac{\underset{d_{\min }}{\overset{d_{\max }}{\int }}n\left(d\right)m\left(d\right)w_{\mathrm{s}}\left(d\right)\mathrm{d}d}{\underset{d_{\min }}{\overset{d_{\max }}{\int }}n\left(d\right)m\left(d\right)\mathrm{d}d}}.\end{array}$$

We only implicitly account for aggregation and fragmentation and explicitly represent the temporally and spatially 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₃, detritus and dust sink in M⁴AGO at the same mean sinking velocity of aggregates (Eq. 3). In contrast to Gehlen et al. (2006), Schwinger et al. (2016) and Heinemann et al. (2019), we explicitly incorporate both variable aggregate size and ballasting through heterogeneous composition. Under the above assumptions, we derive the terms for b, m(d), w_s(d), d_min and d_max in the following sections.

2.2.1 Representation of aggregate microstructure and heterogeneous composition

Marine aggregates are porous (Alldredge and Gotschalk, 1990) 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 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

$$\begin{array}{}\text{(4)}& m\left(d\right)=m_{\mathrm{f}}{d}^{d_{\mathrm{f}}},\end{array}$$

where m_f is a mass factor for the smallest entity. Thus, the density of an aggregate ρ_f decreases with increasing diameter. Aggregates consist of, for example, phytoplankton cells or coccolithophore shells (Alldredge and Gotschalk, 1990), which we consider to be spherical primary particles. Primary particles exhibit their own density ρ_p and diameter d_p. Taking the fractal scaling of aggregate mass into account, the excess density of an aggregate $\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{f}}={\mathit{\rho }}_{\mathrm{f}}-\mathit{\rho }$ with respect to surrounding fluid density ρ can be expressed as (Kranenburg, 1994)

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{\mathit{\rho }}_{\mathrm{f}}=({\mathit{\rho }}_{\mathrm{p}}-\mathit{\rho }){\left({\displaystyle \frac{d_{\mathrm{p}}}{d}}\right)}^{\mathrm{3}-d_{\mathrm{f}}},\text{for }d\ge d_{\mathrm{p}}.\end{array}$$

Furthermore, the aggregate porosity, ϕ, is defined as

$$\begin{array}{}\text{(6)}& \mathit{\varphi }=\mathrm{1}-{\left({\displaystyle \frac{d_{\mathrm{p}}}{d}}\right)}^{\mathrm{3}-d_{\mathrm{f}}},\text{for }d\ge d_{\mathrm{p}},\end{array}$$

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 sinking velocity, w_s:

$$\begin{array}{}\text{(7)}& w_{\mathrm{s}}=\sqrt{{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{f}}}{\mathit{\rho }}}{\displaystyle \frac{gd}{c_{\mathrm{D}}}}}.\end{array}$$

For small particle Reynolds numbers, ${\mathit{\text{Re}}}_{\mathrm{p}}=w_{\mathrm{s}}d/\mathit{\nu }<\mathrm{0.1}$, the drag coefficient is $c_{\mathrm{D}}=\mathrm{24}/{\mathit{\text{Re}}}_{\mathrm{p}}$ and the sinking velocity, w_s, becomes

$$\begin{array}{}\text{(8)}& w_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{\mathrm{18}\mathit{\mu }}}({\mathit{\rho }}_{\mathrm{p}}-\mathit{\rho })g{d}_{\mathrm{p}}^{\mathrm{3}-d_{\mathrm{f}}}{d}^{d_{\mathrm{f}}-\mathrm{1}},\end{array}$$

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 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 and Hill, 2006). With M⁴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, coccoliths, dust particles and detritus as principal components of marine aggregates.

Bushell and Amal (1998) derived a representation of the mean primary-particle size,

$$\begin{array}{}\text{(9)}& \langle d_{\mathrm{p}}\rangle ={\left({\displaystyle \frac{{\sum }_i{n}_i{d}_{\mathrm{p},i}^{\mathrm{3}}}{{\sum }_i{n}_i{d}_{\mathrm{p},i}^{d_{\mathrm{f}}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}-d_{\mathrm{f}}}},\end{array}$$

for an aggregate that is composed of $n_{\mathrm{p}}={\sum }_i{n}_i$ mono-dense spherical primary particles of different diameters d_p,i. Poly-sized 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 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 number of mean, n, and individual primary particles, ∑_in_i (hence, $n\cdot \langle d_{\mathrm{p}}{\rangle }^{\mathrm{3}}={\sum }_i{n}_i{d}_{\mathrm{p},i}^{\mathrm{3}}$ with $n\ne {\sum }_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 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_{\mathrm{p}}=\mathrm{constant}$. This further implies that the ratio K_i, between the total number of primary particles of one particle type, n_i,tot, and the total number of primary particles, $\sum n_{i,\mathrm{tot}}$, in a unit volume is equal to the ratio found in an individual aggregate:

$$\begin{array}{}\text{(10)}& K_i={\displaystyle \frac{n_i}{{\sum }_i{n}_i}}={\displaystyle \frac{n_{i,\mathrm{tot}}}{{\sum }_i{n}_{i,\mathrm{tot}}}}.\end{array}$$

Rewriting Eq. (10), $n_i=K_i{\sum }_i{n}_i$, and inserting it into Eq. (9) gives

$$\begin{array}{}\text{(11)}& \langle d_{\mathrm{p}}\rangle ={\left({\displaystyle \frac{{\sum }_i{K}_i{d}_{\mathrm{p},i}^{\mathrm{3}}}{{\sum }_i{K}_i{d}_{\mathrm{p},i}^{d_{\mathrm{f}}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}-d_{\mathrm{f}}}},\end{array}$$

where the factors K_i can be expressed in HAMOCC via the concentration of each aggregate-forming tracer C_i. Namely, 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-particle diameter, volume $V_{\mathrm{p},i}=\frac{\mathrm{1}}{\mathrm{6}}\mathit{\pi }{d}_{\mathrm{p},i}^{\mathrm{3}}$, and density ρ_p,i:

$$\begin{array}{}\text{(12)}& n_{i,\mathrm{tot}}={\displaystyle \frac{C_i{R}_i}{{\mathit{\rho }}_{\mathrm{p},i}{V}_{\mathrm{p},i}}}.\end{array}$$

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. Ensuring mass conservation, we introduce the volume-weighted primary-particle mean density

$$\begin{array}{}\text{(13)}& \langle {\mathit{\rho }}_{\mathrm{p}}\rangle ={\displaystyle \frac{{\sum }_i{n}_{i,\mathrm{tot}}{V}_{\mathrm{p},i}{\mathit{\rho }}_{\mathrm{p},i}}{{\sum }_i{n}_{i,\mathrm{tot}}{V}_{\mathrm{p},i}}},\end{array}$$

and, hence, multiplication by the volume of the mean primary particle then yields the mass of a mean primary particle (see also Fig. 1):

$$\begin{array}{}\text{(14)}& m(\langle d_{\mathrm{p}}\rangle )={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\mathit{\pi }\langle d_{\mathrm{p}}{\rangle }^{\mathrm{3}}\langle {\mathit{\rho }}_{\mathrm{p}}\rangle .\end{array}$$

Substituting Eq. (14) into Eq. (4), we derive the mass factor for heterogeneous aggregates $m_{\mathrm{f}}=\frac{\mathrm{1}}{\mathrm{6}}\mathit{\pi }\langle d_{\mathrm{p}}{\rangle }^{\mathrm{3}-d_{\mathrm{f}}}\langle {\mathit{\rho }}_{\mathrm{p}}\rangle$. The derivation of the mean primary-particle diameter (Eq. 11), 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_{\mathrm{s}}(d,{\mathit{\rho }}_{\mathrm{p}},d_{\mathrm{p}},\mathrm{\dots })$ can be expressed as $w_{\mathrm{s}}(d,\langle {\mathit{\rho }}_{\mathrm{p}}\rangle ,\langle d_{\mathrm{p}}\rangle ,\mathrm{\dots })$. For a single type of primary particle, all underlying equations reduce to the traditional fractal-scaling relationship (Logan and Wilkinson, 1990; Kranenburg, 1994).

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.

2.2.2 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 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 been suggested as the dominant process in the upper ocean, where turbulent shear reaches typical values of the order of 1 s⁻¹ (Jackson, 1990). By contrast, Alldredge et al. (1990) showed that marine aggregates often withstand oceanic turbulence conditions and suggested biological processes as a 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 of magnitude as ambient turbulent shear in the ocean (Bagster and Tomi, 1974; Adler, 1979; Alldredge et al., 1990). Sinking could thus cause fragmentation in deeper regions of the ocean, where turbulent shear is small (O(0.01 s⁻¹); McCave, 1984; Waterhouse et al., 2014). Since modeling of 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,:

$$\begin{array}{}\text{(15)}& {\mathit{\text{Re}}}_{\mathrm{p}}(d,\langle d_{\mathrm{p}}\rangle ,\langle {\mathit{\rho }}_{\mathrm{p}}\rangle ,d_{\mathrm{f}},\mathit{\nu })={\displaystyle \frac{d{w}_{\mathrm{s}}(d,\langle d_{\mathrm{p}}\rangle ,\langle {\mathit{\rho }}_{\mathrm{p}}\rangle ,d_{\mathrm{f}})}{\mathit{\nu }}}.\end{array}$$

Kiørboe et al. (2001) suggested the particle Reynolds number being in a typical range of up to Re_p=20, while, for example, Alldredge and Gotschalk (1988) measured particle Reynolds numbers of up to Re_p=32. Aggregates thus can exhibit larger Re_p than the laminar case (Re_p<0.1). The drag coefficient, c_D, in Eq. (7) can be represented by the expression for solid spheres of White (2005), valid up to Re_p<10⁵:

$$\begin{array}{}\text{(16)}& c_{\mathrm{D}}={\displaystyle \frac{\mathrm{24}}{{\mathit{\text{Re}}}_{\mathrm{p}}}}+{\displaystyle \frac{\mathrm{6}}{\mathrm{1}+\sqrt{{\mathit{\text{Re}}}_{\mathrm{p}}}}}+\mathrm{0.4}.\end{array}$$

This drag representation leads to smaller settling velocities for large aggregates than the classical Stokes drag ($c_{\mathrm{D}}=\mathrm{24}/{\mathit{\text{Re}}}_{\mathrm{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)

$$\begin{array}{}\text{(17)}& c_{\mathrm{D}}\left({\mathit{\text{Re}}}_{\mathrm{p}}\right)=a_j{\mathit{\text{Re}}}_{\mathrm{p}}^{-b_j},\end{array}$$

to avoid iteration and to allow for an analytical solution of Eq. (3). Applying the parameter values of $a_{j=\mathrm{1}}=\mathrm{24.00},b_{j=\mathrm{1}}=\mathrm{1},\mathrm{for}{\mathit{\text{Re}}}_{\mathrm{p}}\le \mathrm{0.1}$; $a_{j=\mathrm{2}}=\mathrm{29.03},b_{j=\mathrm{2}}=\mathrm{0.871}$, for $\mathrm{0.1}<{\mathit{\text{Re}}}_{\mathrm{p}}\le \mathrm{10}$; and $a_{j=\mathrm{3}}=\mathrm{14.15},b_{j=\mathrm{3}}=\mathrm{0.547}$, for $\mathrm{10}<{\mathit{\text{Re}}}_{\mathrm{p}}\le \mathrm{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 the sinking velocity becomes

$$\begin{array}{}\text{(18)}& w_{\mathrm{s}}={\left({\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}{\displaystyle \frac{\langle {\mathit{\rho }}_{\mathrm{p}}\rangle -\mathit{\rho }}{\mathit{\rho }}}\langle d_{\mathrm{p}}{\rangle }^{\mathrm{3}-d_{\mathrm{f}}}g{\displaystyle \frac{d^{b_j+d_{\mathrm{f}}-\mathrm{2}}}{a_j{\mathit{\nu }}^{b_j}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}-b_j}}.\end{array}$$

By substituting Eq. (18) into Eq. (15), the piecewise integration boundaries, $d_j({\mathit{\text{Re}}}_{\mathrm{p},j=\mathrm{0}\mathrm{\dots }\mathrm{3}}=\mathrm{0},\mathrm{0.1},\mathrm{10},{\mathit{\text{Re}}}_{\mathrm{crit}})$, according to the c_D approximation for Eq. (3), become a function of Re_p:

$$\begin{array}{}\text{(19)}& d_j\left({\mathit{\text{Re}}}_{\mathrm{p},j}\right)={\displaystyle \frac{{\left({\mathit{\text{Re}}}_{\mathrm{p},j}\mathit{\nu }\right)}^{\frac{\mathrm{2}-b_j}{d_{\mathrm{f}}}}}{{\left(\frac{\mathrm{4}}{\mathrm{3}}\frac{\langle {\mathit{\rho }}_{\mathrm{p}}\rangle -\mathit{\rho }}{\mathit{\rho }}\langle d_{\mathrm{p}}{\rangle }^{\mathrm{3}-d_{\mathrm{f}}}g\frac{\mathrm{1}}{a_j{\mathit{\nu }}^{b_j}}\right)}^{\frac{\mathrm{1}}{d_{\mathrm{f}}}}}}.\end{array}$$

Consequently, the concentration-weighted mean sinking velocity (Eq. 3) can then be expressed as

$$\begin{array}{}\text{(20)}& \langle w_{\mathrm{s}}\rangle ={\displaystyle \frac{\sum _{j=\mathrm{0}}^{\mathrm{2}}\left(\underset{max(\langle d_{\mathrm{p}}\rangle ,d_j({\mathit{\text{Re}}}_{\mathrm{p},j}\left)\right)}{\overset{d_{j+\mathrm{1}}\left({\mathit{\text{Re}}}_{\mathrm{p},j+\mathrm{1}}\right)}{\int }}n\left(d\right)m\left(d\right)w_{\mathrm{s}}(d,a_{j+\mathrm{1}},b_{j+\mathrm{1}})\mathrm{d}d\right)}{\underset{\langle d_{\mathrm{p}}\rangle }{\overset{d_{\max }}{\int }}n\left(d\right)m\left(d\right)\mathrm{d}d}},\end{array}$$

where d_max=d_j(Re_crit) is the maximum diameter of aggregates, and by applying $a_{j=\mathrm{0}}=b_{j=\mathrm{0}}=\mathrm{1}$, the lower integration boundary equals the mean primary-particle diameter.

2.2.3 The particle distribution slope, b

Observed aggregate size spectra in the ocean exhibit a spatio-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 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 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 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 and Logan, 1991)

$$\begin{array}{}\text{(21)}& b={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{3}+d_{\mathrm{f}}+{\displaystyle \frac{\mathrm{2}+d_{\mathrm{f}}-min(\mathrm{2},d_{\mathrm{f}})}{\mathrm{2}-b_J}}\right),\end{array}$$

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=\mathrm{2}}$. The assumption of differential-settling-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).

2.2.4 Heuristic approach to variable aggregate 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, 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, 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 aggregation-priming, sticky materials (Alldredge et al., 1993; Azam and 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 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,

$$\begin{array}{}\text{(22)}& \langle \mathit{\alpha }\rangle ={\displaystyle \frac{\mathrm{1}}{A}}\sum _i{n}_i{A}_i{\mathit{\alpha }}_i,\mathrm{where}A=\sum _i{n}_i{A}_i,\end{array}$$

weighted by the primary-particle surfaces $A_i\propto d_{\mathrm{p},i}^{\mathrm{2}}$. We map the mean stickiness to a range between 0 and 1:

$$\begin{array}{}\text{(23)}& \langle \mathit{\alpha }{\rangle }_{\mathrm{map}}={\displaystyle \frac{\langle \mathit{\alpha }\rangle -{\mathit{\alpha }}_{\min }}{{\mathit{\alpha }}_{\max }-{\mathit{\alpha }}_{\min }}},\end{array}$$

where ${\mathit{\alpha }}_{\min }=min\left({\mathit{\alpha }}_i\right)$ and ${\mathit{\alpha }}_{\max }=max\left({\mathit{\alpha }}_i\right)$.

Nicolás-Carlock et al. (2016) introduced a scaling parameter for the effective aggregation range in microscopic aggregation models to stipulate a defined fractal dimension across aggregate sizes. The scaling parameter can be perceived as an indicator of stickiness that defines the effective aggregation range; i.e., higher stickiness results in a larger effective aggregation range. We introduce, as an analogy for the dependency of the fractal dimension on the scaling parameter of Nicolás-Carlock et al. (2016), a transfer function for the mapped mean stickiness to fractal dimension, d_f(〈α〉_map),

$$\begin{array}{}\text{(24)}& d_{\mathrm{f}}(\langle \mathit{\alpha }{\rangle }_{\mathrm{map}})=d_{\mathrm{f},\max }\exp ({\mathit{\beta }}_{\mathrm{f}}\langle \mathit{\alpha }{\rangle }_{\mathrm{map}}),\end{array}$$

with ${\mathit{\beta }}_{\mathrm{f}}=\log (d_{\mathrm{f},\min }/d_{\mathrm{f},\max })$, 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 et al., 1990; Block et al., 1991; Nicolás-Carlock et al., 2016).

2.2.5 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 coccolith. 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), aggregate-incorporated large diatom shells likely enhance the sinking 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 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).

Figure 2 Diatom frustule and the remineralization state-dependent composition of the void. l denotes the thickness of the opal shell with volume V_opal; V_aq and V_POM are the encapsulated volumes of water and POM, respectively. d_p,frustule is the diameter of the diatom frustule.

The opal volume of a modeled diatom is

$$\begin{array}{}\text{(25)}& V_{\mathrm{opal}}={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\mathit{\pi }(d_{\mathrm{p},\mathrm{frustule}}^{\mathrm{3}}-(d_{\mathrm{p},\mathrm{frustule}}-\mathrm{2}l{)}^{\mathrm{3}}),\end{array}$$

where d_p,fustule is the diameter of the diatom, whose opal shell thickness l is expressed in terms of the fixed opal-to-phosphorus formation ratio. The number of diatom frustules per unit volume,

$$\begin{array}{}\text{(26)}& n_{\mathrm{frustule}}={\displaystyle \frac{\left[\mathrm{opal}\right]R_{\mathrm{opal}}}{{\mathit{\rho }}_{\mathrm{opal}}{V}_{\mathrm{opal}}}},\end{array}$$

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 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 hold, it is replaced with the respective volume of water V_aq of density ρ. The frustule density is thus

$$\begin{array}{}\text{(27)}& {\mathit{\rho }}_{\mathrm{frustule}}={\displaystyle \frac{V_{\mathrm{opal}}{\mathit{\rho }}_{\mathrm{opal}}+V_{\mathrm{POM}}{\mathit{\rho }}_{\det }+V_{\mathrm{aq}}\mathit{\rho }}{V\left(d_{\mathrm{p},\mathrm{frustule}}\right)}}.\end{array}$$

During growth and decay, diatoms excrete TEPs which are positively buoyant and possess a density of about ρ_TEP=700 to 840 kg m⁻³ (Azetsu-Scott and Passow, 2004; 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 resolved 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, $m_{\mathrm{potential}}=n_{\mathrm{frustule}}(V_{\mathrm{POM}}+V_{\mathrm{aq}}){\mathit{\rho }}_{\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

$$\begin{array}{}\text{(28)}& {\mathit{\alpha }}_{\mathrm{diatom}}={\displaystyle \frac{m_{\mathrm{e}}}{m_{\mathrm{potential}}}}{\mathit{\alpha }}_{\mathrm{TEP}}+\left(\mathrm{1}-{\displaystyle \frac{m_{\mathrm{e}}}{m_{\mathrm{potential}}}}\right){\mathit{\alpha }}_{\mathrm{opal}},\end{array}$$

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 freshness of detritus. Eventually, the diatom density, ρ_diatom, becomes

$$\begin{array}{}\text{(29)}& {\mathit{\rho }}_{\mathrm{diatom}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{frustule}}{m}_{\mathrm{potential}}+{\mathit{\rho }}_{\mathrm{TEP}}{m}_{\mathrm{e}}}{m_{\mathrm{potential}}+m_{\mathrm{e}}}}.\end{array}$$

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).

2.3 Temperature-dependent opal dissolution and POC remineralization

Marine aggregates tie heterogeneous components together that are disparately remineralized or dissolved. By contrast, in the standard model, detritus, opal and CaCO₃ were sinking separately from each other, and the global remineralization and dissolution rates are tuned independently because the processes are artificially decoupled. In M⁴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 and remineralization.

Opal dissolution is temperature-dependent (Ragueneau et al., 2000, 2006) and is microbially mediated (Bidle et al., 2002). Intact diatom frustules are protected from dissolution by an organic matrix (Lewin, 1961). Once the organic 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 coating of the silicate frustule and (ii) the microbially mediated dissolution of opal with a temperature dependency of Q₁₀≈2.3 (Bidle et al., 2002). We here focus on the temperature-dependent microbially mediated dissolution and introduce a Q₁₀ temperature-dependent opal dissolution:

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\partial }{\partial t}}\left[\mathrm{opal}\right]{|}_{\mathrm{dissolution}}=-r_{\mathrm{opal}}{Q}_{\mathrm{10},\mathrm{opal}}^{\frac{T-T_{\mathrm{ref},\mathrm{opal}}}{\mathrm{10}}}\left[\mathrm{opal}\right],\end{array}$$

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 dissolution (${\partial }_t\left[\mathrm{opal}\right]=-r_{\mathrm{opal}}\left(\mathrm{0.1}\right(T+\mathrm{3}\left)\right)\left[\mathrm{opal}\right]$; Ragueneau et al., 2000; Segschneider and Bendtsen, 2013).

Analogously to opal, we incorporate a temperature-dependent Q₁₀ factor to aerobic POC remineralization (Dell et al., 2011; Mislan et al., 2014) which depends on oxygen concentration, [O₂] (Mauritsen et al., 2019), where $K_{{\mathrm{O}}_{\mathrm{2}}}$ is the half-saturation constant in Michaelis–Menten kinetics, and r_POC is the remineralization rate at reference temperature T_ref,POC:

$$\begin{array}{}\text{(31)}& R_{\mathrm{POC},\mathrm{remin}}=-r_{\mathrm{POC}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}{Q}_{\mathrm{10},\mathrm{POC}}^{\frac{T-T_{\mathrm{ref},\mathrm{POC}}}{\mathrm{10}}}.\end{array}$$

We keep the anaerobic remineralization temperature-independent, 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_{\mathrm{10},\mathrm{POC}}=\mathrm{1}$).

2.4 Model setup, parametrization and evaluation

2.4.1 General model setup

The M⁴AGO sinking scheme was implemented in HAMOCC, which is coupled to the MPIOM (Jungclaus et al., 2013). For the flow of calculations in the M⁴AGO sinking scheme, see Fig. 3. We run both the standard and the M⁴AGO run in a GR15/L40 setup with climatological forcing. This translates to a horizontal resolution of about 1.5^∘ and 40 uneven vertical layers with highest resolution in the first few hundred meters of the ocean. The climatological atmospheric boundary conditions are derived from the second European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis project (ERA-40;  Simmons and Gibson, 2000; Röske, 2005). The mean annual cycle of wind stress, heat and freshwater fluxes is 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₃ due to sedimentation and subsequent burial was accounted for through homogeneously applied weathering rates which were adjusted for the standard run (and the M⁴AGO run): globally, we add ≈99.6 (101.5) G mol P yr⁻¹ as dissolved organic phosphorus and ≈3.2 (2.3) T mol Si yr⁻¹. To compensate for the loss of CaCO₃, we add ≈17.2 (26.5) T mol C yr⁻¹ to surface dissolved inorganic carbon (DIC) and a corresponding amount to surface total alkalinity, A_T, as in DIC : 2A_T. We start the M⁴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 still see drifts of nutrient concentrations in deep, old North Pacific waters at this state (i.e., on average ∼7.1 µmol P m⁻³ 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 throughout the euphotic and mesopelagic zone.

Figure 3 Flow diagram of calculations for the M⁴AGO sinking scheme carried out at every ocean grid point and time step. Marine aggregates in M⁴AGO are composed of spherical primary particles derived from HAMOCC tracers. Primary particles featuring size, density and stickiness are detritus, diatom frustules, coccoliths (CaCO₃) 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.

2.4.2 Parameters of the M⁴AGO scheme

The M⁴AGO sinking scheme introduces a set of new parameters, in particular the primary-particle characteristics, which require constraining and tuning (summarized in Table 1).

Hansen and Kiørboe (1997)Simon et al. (2014)Filella (2007)Kiørboe et al. (1990)Dam and Drapeau (1995)Alldredge and McGillivary (1991)Tambo and Hozumi (1979)Jiang and Logan (1991)Maher et al. (2010)Mahowald et al. (2014)Young and Ziveri (2000)Henderiks (2008)Litchman et al. (2009)Logan and Wilkinson (1990)Kranenburg (1999)Bidle et al. (2002)Mislan et al. (2014)Dell et al. (2011)Bidle et al. (2002)Iversen and Ploug (2010)Azetsu-Scott and Passow (2004)Fettweis (2008)Fettweis (2008)Fettweis (2008)Takahashi et al. (1985)Alldredge and Gotschalk (1988)Kiørboe et al. (2001)

Table 1 Model parameters for M⁴AGO^a and, if adjusted in comparison to the standard HAMOCC, the standard values.

^a Refers to manual tuning of parameter – all other parameters were fixed from beginning. ^b Measured dissolution and remineralization rates include temperature-dependence, which is represented by the Q₁₀ factor in M⁴AGO.

We applied HAMOCC's standard sediment densities for opal, CaCO₃ and dust to the densities of primary particles, namely ρ_opal=2200 kg m⁻³, ρ_calc=2600 kg m⁻³ and ρ_dust=2600 kg m⁻³. For suspended detritus, we chose ρ_det=1100 kg m⁻³ (Fettweis, 2008). As density of TEPs, we applied ρ_TEP=800 kg m⁻³, which is within the measured range of 700 to 840 kg m⁻³ (Azetsu-Scott and Passow, 2004).

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 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 heterogeneity 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 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 (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 $\langle \mathit{\alpha }\rangle \approx O\left({\mathrm{10}}^{-\mathrm{2}}\right)$ (Tambo and Hozumi, 1979). Since stickiness of in situ primary particles is seldom measured, we choose it to our best knowledge and order the stickiness for modeled primary particles according to the mean stickiness of observed aggregate types: ${\mathit{\alpha }}_{\mathrm{dust}}<{\mathit{\alpha }}_{\mathrm{opal}}\le {\mathit{\alpha }}_{{\mathrm{CaCO}}_{\mathrm{3}}}<{\mathit{\alpha }}_{\det }<{\mathit{\alpha }}_{\mathrm{TEP}}$ (see also Table 1). Hence, detritus and TEP-rich aggregates are modeled with a loose structure and low fractal 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 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 (Logan 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₃ are primarily related to phytoplankton mineral structures such 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 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 coccoliths of about 3 to 4 µm in diameter (Young and Ziveri, 2000). We set $d_{\mathrm{p},\mathrm{calc}}=\mathrm{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 (Armbrust, 2009), with a dominant size range of about 12 to 58 µm (equivalent spherical diameter of body volumes 10³–10⁵ µm³; Litchman et al., 2009). We define the diatom frustule size as $d_{\mathrm{p},\mathrm{frustule}}=\mathrm{20}$ µm. The primary-particle size of detritus is weakly constrained and likely ranges from sub-micrometers of microgels and bacteria to millimeter scales of zooplankton body structures (Verdugo et al., 2004). We here chose $d_{\mathrm{p},\det }=\mathrm{4}$ µm. Aeolian dust particles feature a typical 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_{\mathrm{p},\mathrm{dust}}<d_{\mathrm{p},\mathrm{calc}}\le d_{\mathrm{p},\det }<d_{\mathrm{p},\mathrm{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, 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₁₀ temperature-dependent in M⁴AGO. Typically, the Q₁₀ factor for biological processes is in the range between 2 and 3. We here chose $Q_{\mathrm{10},\mathrm{POC}}=\mathrm{2.1}$ (similar to  Mislan et al., 2014, who applied $Q_{\mathrm{10},\mathrm{POC}}=\mathrm{2.0}$). For opal, we tuned $Q_{\mathrm{10},\mathrm{opal}}=\mathrm{2.6}$ compared to 2.3 suggested by Bidle et al. (2002).

2.4.3 Model tuning and evaluation

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 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 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, sediment 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.

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 adjustment 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 primary-particle characteristics, we kept the stickiness values, once 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 remineralization and dissolution rates of POC and opal, respectively. We choose this strategy since reliable remineralization and dissolution 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₃, and their fluxes to sediment within estimated literature ranges. This let us minimally vary the fraction of maximum CaCO₃ production, c_max, in M⁴AGO compared to the standard run.

We here compare and evaluate the M⁴AGO run with respect to (i) where possible, the standard run, (ii) the regional transfer efficiency, as derived by Weber et al. (2016), (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 receive 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.

3 Results and discussion

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⁴AGO, the pattern of POC and associated minerals determines the aggregate properties, which we explicate in 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 biogeochemical tracer distributions (Sect. 3.7). In Sect. 3.8, we discuss the consequence of the transfer efficiency pattern on regional CO₂ 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 assumptions of M⁴AGO (Sect. 3.10).

3.1 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 transport. The global pattern of annual mean integrated primary production therefore remains similar between the standard and the M⁴AGO run (Fig. 4). Globally integrated, the annual net primary production is ≈55.3 Gt C yr⁻¹ in M⁴AGO compared to ≈44.7  Gt C yr⁻¹ in the standard run.

Figure 4 Yearly mean integrated primary production in (a) the standard and (b) the M⁴AGO run. The export efficiency (p ratio) in (c) the standard and (d) the M⁴AGO run.

The ratio between carbon export out of the euphotic zone at depth z₀=100 m and the net primary production, NPP,

$$\begin{array}{}\text{(32)}& p\text{ ratio}\left(z\right)={\displaystyle \frac{\mathrm{POC}\mathrm{flux}\mathrm{at}\mathrm{depth}z=z_{\mathrm{0}}}{\mathrm{NPP}}},\end{array}$$

provides an estimate of how efficient the export is with respect to the net primary production. Globally, about 5.56 and 6.28 Gt C yr⁻¹ is exported out of the euphotic zone in M⁴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⁴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⁴AGO run thus possesses smaller latitudinal variability in the p ratio compared to the standard run.

In comparison to previous estimates on global primary production, both model runs are well within the range of 30 to 70 Gt C yr⁻¹ and show similar pattern of NPP (Carr et al., 2006). The higher NPP in M⁴AGO is due to the enhanced remineralization rates in surface waters, which also lead to 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., 2011, 2012; Buesseler, 1998; Neuer et al., 2002; 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⁴AGO run reduces the bias with respect to previously found low export efficiencies of about 1 % to 10 % (Buesseler, 1998; Neuer et al., 2002).

Figure 5 Yearly mean flux ratios of opal (SiO₂) to POC in (a) the standard and (b) the M⁴AGO run. Yearly mean flux ratios of CaCO₃ to POC in (c) the standard and (d) the M⁴AGO run.

The exported POC is accompanied by biogenic minerals and dust. The export flux ratios between opal and detritus exhibit a clear latitudinal pattern, with high values in the high 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₃-to-detritus flux ratios are confined to equatorial and subtropical regions where silicate depletion favors calcification (Fig. 5c, d). M⁴AGO exhibits a higher CaCO₃-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₃. The spatial distributions of sinking tracers and their ratios prime the marine aggregate properties in M⁴AGO.

3.2 Spatial distribution of marine aggregate properties

The spatial patterns of detritus and mineral fluxes are reflected in the distribution of aggregate properties (cf. pattern 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 density, 〈ρ_p〉, ranges at export depth (100 m) from about 1100 kg m⁻³ in diatom-dominated regions to 1850 kg m⁻³ in the western Pacific, where CaCO₃-to-detritus export ratios are high (Figs. 5b and d and 6a). In the Arctic, some regions harbor 〈ρ_p〉 of maximum 2600 kg m⁻³ where aggregates in our model are dominated by dust particles (Fig. 6a). Particularly in regions of high CaCO₃-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 ranges from about 2 kg m⁻³ in diatom-dominated regions to 35 kg m⁻³ in calcifier-dominated regions (Fig. 6d). Note that we chose Δ〈ρ_f〉_V to account for the increasing porosity 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₃ as ballasting mineral becomes dissolved below the lysocline, and modeled 〈ρ_p〉 decreases again in the deep ocean (Fig. 7a, d). The mean primary-particle size, 〈d_p〉, ranges between the attributed minimum and maximum primary-particle size of tracers, 2 to 20 µm, 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 only little throughout the water column (Fig. 7b), since we assumed invariance of primary-particle size to remineralization, dissolution and other processes.

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⁻¹)) and high (O(>100 m d⁻¹)) sinking velocities. For the mean stickiness, 〈α〉, volume-weighted mean porosity, 〈ϕ〉_V, and the number distribution slope, b, see Fig. A1.

Figure 7 Modeled marine aggregate properties on the Pacific WOA transect P16, which is located at about 150^∘ W. For symbol descriptions, see caption of Fig. 6.

Our simulated d_max is largest in the surface waters of the Southern Ocean and upwelling regions, where TEP-rich aggregates 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 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 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 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 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⁴AGO considers particle sizes ranging from few micrometers to millimeters and thus the full size spectrum, where sinking velocities of O(1 m d⁻¹) to O(>100 m d⁻¹) are represented. 〈w_s〉 thus can significantly differ from reported sinking velocities for large individual aggregates. At the export depth, 〈w_s〉 ranges from about 10 m d⁻¹ in the Southern Ocean and North Pacific region to about 35 m d⁻¹ in the western equatorial Pacific and reaches maximum values of ∼48 m d⁻¹ in dust-dominated Arctic regions (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 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⁻¹. Along the Equator, outside the OMZs, 〈w_s〉 exhibits similarly high values of about 55 m d⁻¹ (Fig. 7f). 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⁻¹) throughout the water column 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⁻¹ (Fig. 7f). Diatom-dominated aggregates in surface waters feature a high buoyancy through TEPs and a loose structure which diminishes with continuous remineralization during their descent.

The aggregate properties entering the M⁴AGO scheme are all directly or indirectly measurable. The comparison of simulated and measured aggregate properties is, however, difficult, as M⁴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 therefore compare M⁴AGO to field and laboratory measurements which examine subsets of the simulated aggregate characteristics.

The modeled aggregate excess densities in diatom-dominated regions compare well to former field and laboratory measurements, where marine aggregates showed excess densities of about 0 to ∼10 kg m⁻³ (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⁻³, which is in the upper range of measured values of 2.1 to 41 kg m⁻³ (Engel et al., 2009). The increased excess density of CaCO₃-dominated aggregates of about 100 to 150 kg m⁻³ in about 1000 m depth compares well with that of observed fecal pellets egested by coccolith-consuming zooplankton (White et al., 2018). These fecal pellets also show similar mean sinking velocities to our modeled 〈w_s〉≈50 m d⁻¹. 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). 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, 1976).

The general difference in typical size between diatom-rich aggregates and CaCO₃ shell-enriched aggregates compares well to observations that also showed smaller CaCO₃-dominated aggregates (Biermann and Engel, 2010). The decreasing 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) 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 balance 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 minerals (Eisma, 1986). In M⁴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 the overall susceptibility of aggregates to shear stress. In M⁴AGO, we disregard this effect and keep Re_crit globally constant.

In summary, the resulting mean sinking velocity in M⁴AGO is of same order of magnitude as that found in 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, 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⁻¹)) as well as rare, but large, fast sinking aggregates (w_s≈O(1000 m d⁻¹)). 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⁴AGO differs significantly from the simple underlying assumption of a linearly increasing sinking velocity with depth in the standard run. M⁴AGO resembles the strongly increasing sinking velocity found in the mesopelagic zone in observations (e.g., Villa-Alfagame et al., 2016) and a modeling approach (DeVries et al., 2014).

3.3 Global pattern of transfer efficiency

The transfer efficiency of POC from z₀=100 m to depth z in the ocean,

$$\begin{array}{}\text{(33)}& T_{\mathrm{eff},\mathrm{POC},z}={\displaystyle \frac{\mathrm{POC}\mathrm{flux}\mathrm{at}\mathrm{depth}z}{\mathrm{POC}\mathrm{flux}\mathrm{at}\mathrm{depth}\mathrm{100}\mathrm{m}}},\mathrm{for}z>z_{\mathrm{0}},\end{array}$$

provides an estimate on the fraction of exported POC that 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,

$$\begin{array}{}\text{(34)}& T_{\mathrm{eff},\mathrm{POC},z}\left(\mathrm{Martin}\right)={\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-\mathit{\beta }},\end{array}$$

to a particular depth and leads to an almost homogeneous global transfer efficiency of $T_{\mathrm{eff},\mathrm{POC},z}\left(\mathrm{Martin}\right)\approx \mathrm{0.12}$ to about 1000 m depth in our standard run (Fig. 8a). The lower remineralization rates in sub-anoxic or even anoxic OMZs compared to oxygenated regions lead to higher transfer efficiencies, visible in the equatorial eastern Pacific Ocean, the northern Indian Ocean, 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 ocean currents and spatially variable turbulent mixing. By contrast, the transfer efficiency in M⁴AGO exhibits a distinct global pattern and possesses higher efficiency in high-latitude and upwelling regions compared to the subtropical gyres where low $T_{\mathrm{eff},\mathrm{POC},z}$ appears. Similar to the standard run, the OMZ regions feature high transfer efficiencies in M⁴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⁴AGO run. We fitted the Martin curve (Eq. 1), to modeled POC fluxes below z₀ to estimate β^′. As expected, the standard run shows little spatial variability apart from the OMZs and features a global, area-weighted mean slope of $\langle {\mathit{\beta }}^{\prime }{\rangle }_{\mathrm{A},\mathrm{Martin}}\approx \mathrm{0.82}$. In agreement with the transfer efficiency pattern, M⁴AGO possesses a strong latitudinal pattern of the effective slope β^′ that varies between ${\mathit{\beta }}^{\prime }\approx$ 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 ${\mathit{\beta }}^{\prime }\approx \mathrm{1.31}$ in the subtropical Pacific gyres (Fig. B1).

Figure 8 Annual mean transfer efficiency for POC in (a) standard and (b) M⁴AGO from export depth (100 m) to about 1000 m (960 m). In (c) the mean transfer efficiency in regions, as defined in (a), is compared to the reconstructed transfer efficiency by Weber et al. (2016). Error bars for Weber et al. (2016) represent uncertainty for the reconstruction of the regional transfer efficiency. For the model results, error bars indicate the spatial standard deviation. For the resulting effective β^′ in M⁴AGO, see Fig. B1. For a seasonal evolution of the transfer efficiency in M⁴AGO, see Fig. C1.

Recently, Weber et al. (2016) reconstructed the global transfer efficiency pattern by diagnosing particulate organic phosphate fluxes. Reconstructing the transfer efficiency via inverse modeling from phosphate concentrations circumnavigates the obstacle of sparse direct observations of fluxes and 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⁴AGO run reveals the inherent inability of the Martin approach to capture the latitudinal variability (Fig. 8a–c). M⁴AGO agrees qualitatively and quantitatively 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⁴AGO run is due to the models overestimation of OMZs extensions (Bopp et al., 2013). The large OMZ 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 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 evolution. 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 M⁴AGO agrees well with the geographical range and pattern found by Buesseler and Boyd (2009) and suggested by Marsay et al. (2015). However, M⁴AGO shows lower maximum values of ${\mathit{\beta }}^{\prime }\approx \mathrm{1.31}$ compared to ${\mathit{\beta }}^{\prime }\approx \mathrm{1.9}$ suggested by Marsay et al. (2015). Globally averaged, M⁴AGO exhibits an effective slope parameter of $\langle {\mathit{\beta }}^{\prime }{\rangle }_{\mathrm{A},{\mathrm{M}}^{\mathrm{4}}\mathrm{AGO}}\approx \mathrm{0.75}$, which is lower than the slope of ${\mathit{\beta }}^{\prime }=\mathrm{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 is in agreement with the elusive trend found by Martin et al. (1987). Overall, M⁴AGO clearly improves the representation of the POC transfer efficiency pattern compared to the standard Martin approach.

3.4 Contributions of 〈w_s〉 and temperature-dependent remineralization to the transfer efficiency pattern

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 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 motions represented by models with this grid resolution. The remineralization length scale (RLS), $z_{\mathrm{POC}}^{*}$, is given by the local ratio of 〈w_s〉 to remineralization (R_POC,remin),

$$\begin{array}{}\text{(35)}& z_{\mathrm{POC}}^{*}={\displaystyle \frac{\langle w_{\mathrm{s}}\rangle }{R_{\mathrm{POC},\mathrm{remin}}}},\end{array}$$

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, β^′.

Figure 9 (a) Remineralization length scale ratio of M⁴AGO to the standard model version at the World Ocean Atlas transect P16. We here focus on showcasing the temperature effect on remineralization and calculated the remineralization length scales, $z_{\mathrm{POC}}^{*}$, without oxygen limitation of remineralization, which cancels out for equal O₂ concentrations. Values smaller than 1 imply stronger POC flux attenuation in M⁴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_{\mathrm{s},\mathrm{POC}}^{*}$, with depth (shown bottom left). Contour lines provide the Q₁₀ factor temperature-dependent remineralization rates with r_POC at reference temperature T_ref in M⁴AGO. In the standard run, the aerobic rate is globally constant (0.026 d⁻¹). (b) Relative contributions of sinking, ${\mathrm{RC}}_{\langle w_{\mathrm{s}}\rangle }$, and remineralization, RC_remin, to difference between standard and M⁴AGO run. Contour lines provide the relative contribution of remineralization.

The RLSs in surface waters and the upper mesopelagic zone of subtropical and equatorial regions are shorter in the M⁴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⁴AGO compared to the standard run in these regions (Fig. 4b). In the mesopelagic zone, the RLSs in M⁴AGO are 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⁴AGO compared to the standard run. In order to analyze which of the 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⁴AGO and the standard run as

$$\begin{array}{}\text{(36)}& {\mathrm{RC}}_{P_i,\mathit{x}}\left[\mathit{\%}\right]={\displaystyle \frac{{\partial }_{P_i}{z}_{\mathrm{POC}}^{*}\cdot \mathrm{\Delta }{P}_{i,\mathit{x}}}{{\sum }_i|{\partial }_{P_i}{z}_{\mathrm{POC}}^{*}\cdot \mathrm{\Delta }{P}_{i,\mathit{x}}|}}\cdot \mathrm{100},\end{array}$$

where the processes, $P_i=\mathit{\{}\langle w_{\mathrm{s}}\rangle ,R_{\mathrm{POC},\mathrm{remin}}\mathit{\}}$, refer to Eq. (35), ${\partial }_{P_i}$ is the partial derivative of $z_{\mathrm{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⁴AGO run and the standard run at spatial point x. M⁴AGO possesses generally higher remineralization rates than the standard run, which would increase the flux attenuation compared to the standard run. The temperature-dependent remineralization in M⁴AGO shows lower rates in the cold waters of the high latitudes than in the equatorial and subtropical regions (Fig. 9a). Within M⁴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 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₃-dominated subtropical gyres, 〈w_s〉 falls below the sinking velocity of the standard model in regions deeper than 3000 m and thus contributes further to the stronger flux attenuation compared to the standard run. The higher 〈w_s〉 in M⁴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 remineralization in M⁴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.

3.5 Impact of mineral size and ballasting effect on sinking velocity

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 understand the underlying factors that control the mean sinking velocity. We suggest in the following that the size of primary particles might be as important as the density of the ballasting material for defining the sinking velocity of aggregates.

M⁴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:

$$\begin{array}{}\text{(37)}& R_{X_i,z}\left[\mathrm{\%}\right]={\displaystyle \frac{{\partial }_{X_i}\langle w_{\mathrm{s}}\left({\stackrel{\mathrm{‾}}{X}}_{i,z}\right)\rangle \cdot \mathrm{\Delta }{X}_{i,z}}{{\sum }_i\left|{\partial }_{X_i}\langle w_{\mathrm{s}}\left({\stackrel{\mathrm{‾}}{X}}_{i,z}\right)\rangle \cdot \mathrm{\Delta }{X}_{i,z}\right|}}\cdot \mathrm{100},\end{array}$$

where $X_i=\mathit{\{}\langle {\mathit{\rho }}_{\mathrm{p}}\rangle ,\langle d_{\mathrm{p}}\rangle ,\mathit{\mu },d_{\max },b,d_{\mathrm{f}}\mathit{\}}$, $\mathrm{\Delta }{X}_{i,z}=X_{i,z}-{\stackrel{\mathrm{‾}}{X}}_{i,z}$ and $\sum R_{X_i,z}\ne \mathrm{100}$, but $\sum \left|R_{X_i,z}\right|=\mathrm{100}$, and ${\stackrel{\mathrm{‾}}{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 percentage-wise contribution to the local sinking velocity by local primary-particle density compared to the global mean primary-particle density? By neglecting the higher-order terms, we provide only qualitative insights into the role of the different aggregate properties and molecular viscosity in 〈w_s〉.

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.

CaCO₃ 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 contributes positively to sinking velocity (Fig. 10d). Additionally, in our model, d_f increases particle excess density and thus 〈w_s〉 in CaCO₃-rich areas (Fig. 10c). In the Southern Ocean, North Pacific, North Atlantic and upwelling regions, particularly the large 〈d_p〉 in diatom-dominated regions enhances 〈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 exogenous 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 contrasting relative contributions of dense mean primary particles and large mean primary particles to 〈w_s〉 become even more pronounced. Denser, small CaCO₃ particles contribute to 〈w_s〉 in the subtropical, oligotrophic regions of the oceans, while comparably less dense but larger opal frustules lead to high 〈w_s〉 in nutrient-rich upwelling and high-latitude regions. Even though our formulation for 〈w_s〉 is highly non-linear, 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 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 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.

Microstructure d_f directly prescribes the aggregate number distribution slope, b, in M⁴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⁴AGO, where b varies between ≈ 3.19 and 3.76 and thus likely underestimates the spatial variability and relative contribution of b to 〈w_s〉. At present, M⁴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 required 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₃ acts as a strong ballasting agent (e.g., Francois et al., 2002; Balch 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 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 sinking 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 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, 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., 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 morphology effects on POC fluxes. The higher size and morphology variability in diatom-dominated phytoplankton communities, 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₃ is weaker 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.

3.6 Regional fluxes and rain ratios

Biogenic minerals, dust particles and detritus are tied together in M⁴AGO and define the composition, microstructure and thus sinking velocity of aggregates. In turn, the tracers are independently remineralized or dissolved. Both processes 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 opal (and CaCO₃ 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, ${\mathit{\beta }}_{\mathrm{opal}}^{\prime }$, exhibits similar regional variability in the M⁴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 ${\mathit{\beta }}_{\mathrm{opal}}^{\prime }$. Similar to POC fluxes, 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 velocity than the 25 m d⁻¹ 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).

Figure 11 Area-weighted mean effective power function slopes for opal fluxes, $\langle {\mathit{\beta }}_{\mathrm{opal}}^{\prime }{\rangle }_{\mathrm{A}}$, in the regions defined in Fig. 8. Calculated for opal fluxes below 100 m. Error bars show regional standard deviation.

In sum, M⁴AGO estimates that globally ∼1.03 Gt Si yr⁻¹ (∼1.04 Gt Si yr⁻¹ in the standard run) reaches the seafloor, which is in good agreement with present estimates of about 1 Gt Si yr⁻¹ (Tréguer, 2002). Generally, the modeled ${\mathit{\beta }}_{\mathrm{opal}}^{\prime }$ is well within the range of observations with estimates of ${\mathit{\beta }}_{\mathrm{opal}}^{\prime }=\mathrm{0.22}\pm \mathrm{0.53}$ (Boyd et al., 2017). In M⁴AGO, dissolution of opal has the implication that ballasting mineral ratios between opal and CaCO₃ shift towards higher importance of CaCO₃ with depth. Since the dissolution of opal is temperature-dependent, the preservation efficiency for opal-to-CaCO₃ 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 therefore 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 ranging 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 at their respective depths. The M⁴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⁴AGO compared to the standard run (Fig. 12).

Figure 12 Monthly POC∕Si rain ratios in the standard run and in M⁴AGO compared to the monthly climatological mean derived from the Mouw et al. (2016a, b) data set. The black line denotes the 1 : 1 line. Notice that the axes are in log and have different limits among regions but are comparable between runs. No data are available for opal fluxes in the equatorial tropical Atlantic (ETA), where opal fluxes are small (refer to Fig. 5a, b). Correlation coefficient, r, and significance value, p, are given if r>0.4.

However, M⁴AGO introduces a bias in regions deeper than ∼4000 m towards smaller POC∕Si ratios. The RLSs in the deep ocean are hence too short in M⁴AGO. While the scatter of the POC∕Si ratio reduces in M⁴AGO, the overall variability in the flux ratios becomes compressed compared to the measured variability. This compression might be caused by 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 silicification 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 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., 2003, and references therein).

3.7 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 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, 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 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 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 in the mesopelagic zone, 362 m.

Figure 13 (a–d) Taylor diagram for tracers (oxygen – O₂; nitrate – ${\mathrm{NO}}_{\mathrm{3}}^{-}$; phosphate – ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$; 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). (a, c) Pacific and (b, d) Atlantic for (a, b) 6 and 100 m; (c, d) are for 362 and 960 m. (e) Surface phosphate concentration in WOA (f) in the standard run and (g) difference plot between M⁴AGO and the standard run.

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 distributions 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 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 contribute 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 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 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 latter tracers are therefore generally more prone to biogeochemical model reformulation, and larger differences can be expected between M⁴AGO and the standard run. In addition, M⁴AGO links the cycles of phosphate and silicate closer to nitrate and oxygen through common sinking of particulate matter. Phosphate and silicate therefore experience an additional biogeochemical cycle imprint in M⁴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 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. 13e–g).

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 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 currents, 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⁴AGO, the water column above 800 m is more oxygenated and below 800 m 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 manner; i.e., denitrification and sulfate reduction increase alkalinity, whereas aerobic remineralization decreases alkalinity. While A_T and silicate are generally worse than other tracers when compared to observations, M⁴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 compared 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⁴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 also increases the retention time of the tightly coupled opal in M⁴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⁴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⁴AGO, since the silicate distribution in the euphotic zone directly affects opal production and thus the sinking velocity of aggregates, transfer efficiency and nutrient distributions. The tighter coupling of silicate to other biogeochemical cycles in M⁴AGO therefore warrants further future investigation.

3.8 Regional CO₂ fluxes

In a 100-year climatological steady state, both model runs show the Southern Hemisphere ocean acting as a net source of CO₂ to the atmosphere, while the Northern Hemisphere acts as a net sink (Fig. 14a). Consequently, net oceanic CO₂ transport across the Equator from the Northern to the Southern Hemisphere exists. In the simulation with M⁴AGO, a stronger CO₂ 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₂ to the atmosphere. As is visible from the decreasing difference between the cumulative zonal CO₂ fluxes, M⁴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⁻¹ between the two model runs remains, the small residual CO₂ flux imbalance of 0.07 and 0.05 Gt C yr⁻¹ in M⁴AGO and the standard run, respectively, indicate well-spun-up model runs. In general, the latitudinal zonal CO₂ fluxes and the cross-equatorial southward oceanic CO₂ transport agree qualitatively and quantitatively well with former forward-integrated models for preindustrial conditions (e.g., Sarmiento et al., 2000; Gloor et al., 2003; Mikaloff Fletcher et al., 2007).

Figure 14 (a) Climatological cumulative zonal CO₂ flux in the standard and the M⁴AGO run (from south to north). Generally, negative fluxes represent net CO₂ uptake by the ocean. Note the second y axis for the difference between the two runs (also in c–k). (b) Global annual net CO₂ flux under linearly increasing atmospheric CO₂ concentration. (c–k) Regional cumulative CO₂ fluxes under increasing atmospheric CO₂ (regions as defined in Fig. 8).

To study the effect of the changed transfer efficiency in M⁴AGO on the CO₂ uptake without indirect climate effects, we perform model runs in which we linearly increase the atmospheric CO₂ 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, ocean temperature, stratification, and thus primary production and remineralization occurs in these model runs.

During the atmospheric CO₂ increase, both model runs show similar global annual oceanic CO₂ uptake (Fig. 14b). About 181–185 Gt C is taken up by the global ocean throughout the 150-year period. The discrepancy of ≈4 Gt C between the two runs arises to a large extent from the residual imbalance of 0.02 Gt C yr⁻¹ uptake, cumulated over 150 years.

The general trends in the regional cumulative fluxes remain the same between the two model runs (Fig. 14c–k), and most regions act as net sinks of CO₂. Only the upwelling regions of the equatorial tropical Atlantic (ETA) and Pacific (ETP) are net sources of CO₂ to the atmosphere at the end of the 150-year period. The magnitude of the regional cumulative fluxes, however, varies among the two model runs.

Since we can rule out temperature-driven effects on the differences in CO₂ uptake when comparing the two climatologically forced runs, a different state and/or changes in A_T and dissolved inorganic carbon (DIC) concentration determine the partial pressure of CO₂, pCO₂ and thus CO₂ fluxes. In our runs, major changes in A_T and DIC can locally occur due to a change in, first, the CaCO₃-to-detritus rain ratio as a result of varying CaCO₃ and POM production and remineralization, or, second, an increase in the POC transfer efficiency, for example, which reflects a transfer of carbon to deeper waters where CO₂ 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₂ 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₂ is lower in M⁴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 similar CO₂ 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₂ fluxes (e.g., in the SAZ). In the STP region, where high CaCO₃-to-detritus rain ratios occur (refer to Fig. 5), A_T is about 0.2 % lower in M⁴AGO than in the standard run, which enhances pCO₂ by about 2 % and thus decreases the CO₂ uptake.

In summary, under linearly increasing CO₂ in a non-interactive mode, the global CO₂ fluxes remain similar for the standard and the M⁴AGO run. However, regional CO₂ fluxes change, which is potentially linked to the changing pattern of transfer efficiency and the CaCO₃-to-detritus rain ratio. A detailed study on the feedbacks between the transfer efficiency, represented by M⁴AGO, and a transient climate on oceanic CO₂ uptake is out of the scope of this paper and will be part of future investigations.

3.9 Sensitivity analysis

With M⁴AGO, the processes of temperature-dependent remineralization and a number of new model parameters to represent 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 efficiency, (ii) phosphate as an essential nutrient for primary production, and (iii) silicate as a nutrient and circulation-reflecting agent.

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⁻¹), 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.

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 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_{\mathrm{p},\mathrm{frustule}}=$20 µm to $d_{\mathrm{p},\mathrm{frustule}}=$10 µm. In the second sensitivity experiment, S(d_f), the varying fractal dimension, ranging between 1.6 and 2.4 in the M⁴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 tuning the M⁴AGO run. Here, we increase the half-saturation constant for phosphate uptake by diazotrophs from 0.05 to 0.10 µmol L⁻¹. 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 the M⁴AGO run reduces 〈w_s〉 and thus the RLS in diatom-dominated regions. As a consequence, the transfer efficiency in silicifier-dominated regions is lower in S(d_p,frustule) than in M⁴AGO (Fig. 15g). Accordingly, opal dissolves closer to surface waters and the silicate concentration increases compared to the M⁴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 increases 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 producers 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⁴AGO, changes of the phytoplankton size structure with ongoing ocean warming could affect RLSs, transfer efficiency 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 references therein). S(d_p,frustule) shows that decreasing primary-particle 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₂ uptake under ongoing 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 calcifier-dominated waters compared to the M⁴AGO run (cf. Fig. 6c). In the M⁴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 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 decreased 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 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⁴AGO (Fig. 15b). As a consequence of the potential sensitivity 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 nitrate on phosphate. They are regionally confined to warm tropical and subtropical regions (Paulsen et al., 2017, 2018). Diazotrophs modulate the phosphate concentration in the subtropical gyres. The diazotrophs' global primary production increases from about 2.2 Gt C yr⁻¹ in the standard run to ≈2.9 Gt C yr⁻¹ in the M⁴AGO run. A slower growth response to phosphate concentrations, enforced by raising the half-saturation constant from 0.05 to 0.10 µmol L⁻¹ in S(K_P,diaz), increases the phosphate concentrations in the subtropical gyres to more than 150 % (Fig. 15c). The primary production through diazotrophs reduces significantly from a former global value of ≈2.9 to ≈1.9 Gt C yr⁻¹ 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⁻¹ in S(d_p,frustule) and S(d_f), respectively, 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 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 to the standard run, this phenomenon is generally intensified in M⁴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⁴AGO, require further future evaluation.

3.10 Current limitations of M⁴AGO

Developing M⁴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 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⁴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 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 fragmentation 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 b≈2 to 5 (Guidi et al., 2009). As a consequence, M⁴AGO probably underestimates the variability in sinking velocity, particularly in the upper ocean, where the size distribution dynamically evolves.

In M⁴AGO, TEPs are simplistically considered. In agreement with Mari et al. (2017), TEPs are microstructure-loosening and buoyancy-adding for diatom-dominated aggregates. For CaCO₃, we only considered the coccolith size which decreases sinking velocity compared to an aggregate composed of intact coccospheres. For simplicity, we 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 fragmentation 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 aggregation (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 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).

4 Summary and conclusions

We implemented the novel Microstructure, Multiscale, Mechanistic, Marine Aggregates in the Global Ocean (M⁴AGO) scheme in HAMOCC to improve the representation of the biological carbon pump in an Earth system model framework. M⁴AGO accounts for the heterogeneity and microstructure of aggregates and thus clearly defines measurable statistic aggregate properties in HAMOCC. M⁴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₁₀ temperature-dependent dissolution and aerobic remineralization of opal and POC, respectively.

In contrast to the standard HAMOCC version, M⁴AGO provides a mechanistic understanding for the recently published global transfer efficiency pattern and represents it 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, co-determine the high POC transfer efficiency in high latitudes and upwelling regions. Our model results support previous findings that CaCO₃ 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₃-to-detritus rain ratio alters regional CO₂ fluxes, while the global uptake remains the same as in the standard run, when atmospheric CO₂ is linearly increased without climate feedbacks.

The highest uncertainties in parametrizing M⁴AGO are with respect to the weakly constrained primary-particle surface 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 development for the representation of marine aggregates would highly benefit from sub-aggregate-scale measurements of microstructure, adhesive surface properties and primary-particle composition.

Our findings and the underlying model concept suggest 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 extended 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 structures in POC fluxes, particularly in nutrient-rich upwelling regions, where a wide, variable size spectrum of diatoms prevails. Second, the indirect temperature effect on phytoplankton 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₂ in the atmosphere.

In conclusion, M⁴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.

Appendix A: Additional aggregate properties

In addition to the aggregate properties presented in Fig. 6, M⁴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 aggregates (Eq. 6) can be deduced from aggregate size, primary-particle 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 therefore 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, 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⁴AGO (Fig. A1). At the bottom of the mesopelagic zone, 〈α〉 is almost homogeneous apart from the OMZ regions, where the lower 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 decline in calcifier-dominated surface waters (Fig. A1c, f). In M⁴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-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 occurrence in deep OMZs and thus their influence on aggregate properties.

Figure A1 Mean stickiness 〈α〉 (Eq. 22), volume-weighted porosity, 〈ϕ〉_V, and the number distribution slope, b (Eq. 21).

Appendix B: Effective Martin slope in M⁴AGO

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 et al., 2015). In M⁴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 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, we found, apart from OMZs, a rather homogeneous β^′ with a global value of $\langle {\mathit{\beta }}^{\prime }{\rangle }_{\mathrm{A},\mathrm{Martin}}\approx \mathrm{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 $\langle {\mathit{\beta }}^{\prime }{\rangle }_{\mathrm{A},\mathrm{Martin}}$ is reduced by the lower anaerobic 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 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.

Figure B1 Effective Martin curve slope for POC in M⁴AGO, β^′.

Appendix C: Seasonal transfer efficiency

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 variation, and so does transfer efficiency (Lutz et al., 2007). For example, if we assume an average sinking speed of about 25 m d⁻¹, 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). Accordingly, 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 transfer efficiency can thus heavily deviate from the climatological state. This seasonal behavior is reflected in M⁴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 similar qualitative pattern, the seasonal amplitude of the transfer efficiency in high latitudes is lower in M⁴AGO compared to the standard run.

Figure C1 Seasonal evolution of the transfer efficiency in M⁴AGO. Note the difference of the color bar values compared to Fig. 8.

Appendix D: Mathematical symbols

Table D1 Mathematical symbols and their description.