The effect of marine aggregate parameterisations on nutrients and oxygen minimum zones in a global biogeochemical model

Particle aggregation determines the particle flux length scale and affects the marine oxygen concentration and thus the volume of oxygen minimum zones (OMZs) that are of special relevance for ocean nutrient cycles and marine ecosystems, and that have been found to expand faster than can be explained by current state-of-the-art models. To investigate the impact of particle aggregation on global model performance, we carried out a sensitivity study with different parameterisations of marine aggregates and two different model resolutions. Model performance was investigated with 10 respect to global nutrient and oxygen concentrations, as well as extent and location of OMZs. Results show that including an aggregation model improves the representation of OMZs. Moreover, we found that besides a fine spatial resolution of the model grid, the consideration of porous particles, an intermediate to high particle sinking speed and a moderate to high stickiness improve the model fit to both, global distributions of dissolved inorganic tracers and regional patterns of OMZs, compared to a model without aggregation. Our model results therefore suggest that improvements not only in the model 15 physics, but also in the description of particle aggregation processes can play a substantial role in improving the representation of dissolved inorganic tracers and OMZs on a global scale. However, dissolved inorganic tracers are apparently not sufficient for a global model calibration, which could necessitate global model calibration against a global observational dataset of marine organic particles.

One potential parameter affecting distributions of dissolved oxygen and thereby the volume and location of OMZs is the biological carbon pump (Volk and Hoffert, 1985). Global ocean model studies show that the biological pump is important for the distribution of dissolved inorganic tracers in the ocean Primeau, 2006, 2008) as well as atmospheric pCO 2 (Kwon et al., 2009;Roth et al., 2014). It further affects the feeding of deep sea organisms (Kiko et al., 2017) as well as the OMZs volume (Kriest and Oschlies, 2015). The biological carbon pump can be subdivided into three components: 5 production of organic matter and biominerals in the euphotic surface layer, particle export into the ocean interior and finally their decomposition in the water column and on the sea floor (Le Moigne et al., 2013). Estimates of the export of organic carbon out of the surface layer range from 5 to 20 Gt C yr -1 , with the large uncertainty illustrating the gap in our understanding of this process (Henson et al., 2011;Honjo et al., 2008;Keller et al., 2012;Laws et al., 2000;Oschlies, 2001).
Further uncertainties are associated with the exact shape of the particle flux profile (e.g. exponential function vs. power law; 10 Banse, 1990;Berelson, 2002;Boyd and Trull, 2007;Buesseler et al., 2007;Lutz et al., 2002;Martin et al., 1987) and its possible variations in space and time. Recent studies suggest conflicting evidence with regard to the spatial variation of the particle flux length scale (Guidi et al., 2015;Marsay et al., 2015), which may again be influenced by the methodology of estimating the particle flux profile and thus the potential sensitivity to the considered depth (Marsay et al., 2015). Also, the underlying mechanisms for a potential spatio-temporal variation remain unclear: some studies attribute this to variations in 15 temperature and associated temperature-dependent variation in remineralisation (Marsay et al., 2015), while other studies derive this from variations in particle size distributions (Guidi et al., 2015).
One mechanism, that leads to a variation in particle size distribution, consists in the formation of marine aggregates, which exhibit variable sinking speeds. For example, Alldredge and Gotschalk (1988) and Nowald et al. (2009) found sinking rates for aggregates ranging between 10 and 386 m d -1 . Particle sinking speed, and thus the particle flux profile, depends on 20 mineral ballast (Armstrong et al., 2002;Ploug et al., 2008), porosity and particle size (Alldredge and Gotschalk, 1988;Kriest, 2002;Smayda, 1970). Large particles are associated with high sinking speed and fast passage through the water column, resulting in low remineralisation and thus a small OMZs volume, and vice versa. It can therefore be expected that particle aggregation favouring fast sinking speeds can alter the volume of OMZs compared to small particles with low sinking speeds (Kriest and Oschlies, 2015). 25 However, there are still some gaps in our understanding of the parameters that control the aggregation rate as well as the particle's sinking behaviour. For example, in-situ measurements show almost no dependency between diameter and sinking speed (Alldredge and Gotschalk, 1988), whereas aggregates produced on a roller table show a noticeable relationship (Engel and Schartau, 1999). Furthermore, values for stickiness, which defines the probability that after collision two particles stick together, vary over a wide range. Stickiness depends on the chemistry of the particle's surface (Metcalfe et al., 2006) and the 30 particle type (e.g. Hansen and Kiørboe, 1997) and ranges between almost zero and one (e.g. Alldredge and McGillivary, 1991;Kiørboe et al., 1990). Thus, aggregation as one process that induces variations in particle size, and thus sinking speed, is only loosely constrained through its parameters.
To explore these relationships further and to examine whether a spatially variable sinking speed improves the fit of a global biogeochemical model to global distributions of dissolved inorganic tracers and regional patterns of OMZs, this study uses the three-dimensional Model of Oceanic Pelagic Stoichiometry (Kriest and Oschlies, 2015), coupled with a module for particle aggregation and size-dependent sinking (Kriest, 2002). Given the large uncertainty associated with parameterisations of marine aggregates, we carried out 36 sensitivity experiments, in which we varied parameters relevant for particle aggregation and sinking. As in previous studies, the model`s fitness is evaluated by the Root Mean Square Error (RMSE) 5 against observational data of dissolved inorganic tracers, namely PO 4 , NO 3 and O 2 . This study additionally determines the model fitness with respect to extent and location of OMZs, following the approach by Cabré et al. (2015).
To examine the above-mentioned questions, and explore the effects and uncertainties of a model that simulates particle dynamics on a global scale for a seasonally cycling stationary ocean circulation, our main questions are as follows: 10 1. Does a model that includes explicit particle dynamics improve the representation of observed PO 4 , NO 3 and O 2 ?
2. Does a model that includes explicit particle dynamics improve the representation of observed OMZs, and do the 'best' parameters with respect to this metric agree with those constrained by dissolved inorganic tracers?
3. What are the effects of uncertainties in the parameterisation of organic aggregates on model results ? 15 4. Can the assumptions inherent in the model confirm either of the spatial particle flux length scale maps proposed by Marsay et al. (2015) or Henson et al. (2015) and Guidi et al. (2015)?
This paper is organised as follows: we first describe the model and its assessment with regard to dissolved inorganic tracers and OMZs, including the sensitivity experiments carried out with the model. We then present the outcome of the sensitivity 20 experiments, with special focus on the metrics defined above. We finally examine and discuss derived maps of particle flux length scales against the background of maps derived from observed quantities Marsay et al., 2015;Guidi et al., 2015).

Oceanic transport 25
In this study, we used the 'Transport Matrix Method' (TMM) (Khatiwala et al., 2005), as an efficient offline method to simulate biogeochemical tracer transport with monthly mean transport matrices (TMs). Additional fields of monthly mean wind, temperature and salinity extracted from the underlying circulation model are used to simulate air-sea gas exchange of oxygen, and to parameterise temperature-dependent growth of phytoplankton. For our experiments, we used two different types of TMs and forcing fields: one set derived from a coarser resolution (hereafter called MIT2.8), and one from a finer 30 resolution version, based on a data-assimilated circulation (ECCO1.0) (Stammer et al., 2004). The MIT2.8 forcing and transport represent a resolution of 2.8° x 2.8° and 15 depth layers with a thickness ranging between 50 m and 690 m. ECCO1.0 TMs and forcing are based on a resolution of 1° x 1° and 23 depth layers, with a thickness ranging between 10 m and 500 m. Further details about the two setups can be found in Kriest and Oschlies (2013).
In general, we used a time step length of 1/2 day for physical transport, and a time step length of 1/16 day for biogeochemical interactions in the coarse resolution, MIT2.8. Because some parameter configurations allow a very large 5 particle sinking speed, which may exceed more than one box per time step, in MIT2.8 we used a biogeochemical time step length of 1/70 day for all simulations with = 1.17 (see Table 1), in the finer resolution, ECCO1.0, we used in all experiments a time step of 1/80 day (see Table 1) but with exception of three experiments, where we used a length of 1/160 day (these are the experiments for a strong increase of sinking speed with particle size, given by parameter = 1.17; see Table 1). Each model was integrated for 3,000 years until tracers approached steady state. The last year is used for analysis 10 as well as misfit calculations.

Model of Oceanic Pelagic Stoichiometry
The Model of Oceanic Pelagic Stoichiometry, called MOPS (Kriest and Oschlies, 2015), is based on phosphorus, and simulates phosphate, phytoplankton, zooplankton, dissolved organic phosphorus (DOP) and detritus. The unit of each tracer 15 is given in mmol P m -3 . In addition, MOPS simulates oxygen and nitrate. The P-cycle is coupled to oxygen by using a fixed stoichiometry of R -O2:P =171.739, and to nitrogen by R P:N =16.
The stoichiometry of anaerobic and aerobic remineralisation is parameterised following Paulmier et al. (2009).
Remineralisation of detritus and DOM is fixed to a constant nominal remineralisation rate r and is dependent on oxygen but independent of temperature. If oxygen concentrations decrease, denitrification replaces aerobic respiration, consuming 20 nitrate. If neither oxygen nor nitrate is sufficiently available, remineralisation stops as the model does not account for other electron acceptors such as sulphate. As both forms of remineralisation follow a saturation curve (Monod-type), the realised remineralisation rate may diverge from the constant nominal remineralisation rate.
On long timescales, the loss of fixed nitrogen through denitrification is balanced by temperature-dependent nitrogen fixation. Therefore, it should be noted that while phosphorus is conserved, the inventory of fixed nitrogen as well as oxygen is 25 variable, and dependent on ocean circulation and biogeochemistry (Kriest and Oschlies, 2015).
In the basic model without aggregation the sinking speed of detritus increases linearly with depth. With constant remineralisation rate r, the particle flux can thus be described by ∝ %& with = ) * (Kriest and Oschlies, 2008), and is therefore (for constant r, e.g. in a fully oxic water column) comparable to the common power-law description of observed particles fluxes (Martin et al., 1987). The fraction of detritus reaching the seafloor follows two pathways: One fraction is re-30 suspended back into the deepest box of the water column and the other one is buried into the sediment and therefore responsible for P-removal. However, the P-budget remains annually unchanged by the resupply of buried P via river runoff.

Model for particle aggregation and size dependent sinking
Different approaches have been applied to simulate particle aggregation in the marine environment. A detailed representation of the particle size spectrum can be accomplished by explicitly simulating many different size classes, which interact with each other via collision-based aggregation, particle sinking, remineralisation and breakup (Burd, 2013;Jackson, 1990). This flexible approach captures the details of the size spectrum and its spatio-temporal variation in a very detailed way. However, 5 it is computationally expensive, and thus prohibitive to be applied to large spatial and long temporal scales.
The aggregation module applied in MOPS parameterises a continuous log-log-linear size distribution of particles via the spectral slope e calculated from number and mass of particles (Kriest and Evans, 2000). The particle size distribution is influenced by size-dependent particle aggregation and sinking (Kriest, 2002;Kriest and Evans, 2000). Because aggregation reduces particle numbers (but not mass), and sinking preferentially removes large particles, number and mass change 10 independently. By assuming a log-log-linear size spectrum, the slope e of this spectrum can, at each time step and grid point, be computed from the particle number and total particle mass.
The model requires parameters for the power-law relationships between particle diameter, d, and mass, m, (m = Cd ζ ) and between particle diameter and sinking speed, w, (w = Bd ) to be specified. In our model experiments, we assign fixed values for the minimum diameter and mass of a primary particle of size of d 1 =0.002 cm and m 1 =0.00075 nmol P. The exponent for 15 the relationship between size and mass is set to ζ = 1.62, as proposed for marine aggregates in Kriest (2002), which is in line with more recent findings (Burd and Jackson, 2009;Jouandet et al., 2014). For the relationship between size and sinking speed we test two alternative values for eta, namely = 0.62 and = 1.17 for the exponent, and w 1 between 0.7-2.8 m d -1 for the minimum sinking speed (see below). Assuming a constant degradation rate, the average sinking speed of all particles combined would increase with depth due to higher sinking speed of large particles and their higher proportion in the deeper 20 ocean interior. To prevent instabilities at very large sinking speeds (very flat size distributions), as in Kriest and Evans (2000) and Kriest (2002) we restrict the size dependency of sinking and aggregation to a maximum diameter of D L . Beyond D L , these processes do not vary with particle size any more. In our model experiments, we let this parameter vary between 1, 2 and 4 cm.
Changes in the number of marine particles are dependent on particle aggregation, described by the collision rate, and the 25 probability that two particles stick together, a. In our model experiments we vary a between 0.2-0.8. The collision rate depends on turbulent shear and differential sinking and is parameterised as in Kriest (2002). We assume that the turbulent shear is high in the euphotic layers and zero in the deeper ocean layers.
To avoid complications and non-linear feedbacks, in the experiments presented here, we assume that plankton mortality and zooplankton egestion as well as quadratic zooplankton mortality produce new detritus particles, but do not change the size 30 spectrum.
By using this setup, the module is similar to parameterisations of particle size applied in other large-scale or global models (Gehlen et al., 2006;Oschlies and Kähler, 2004;Schwinger et al., 2016).

Adjustment of biogeochemical model parameters 5
Introducing aggregates and a dynamic particle flux profile to the global model MOPS has a strong impact on biogeochemical model dynamics. Starting from parameter values of the calibrated model setup (without aggregation) of , we calibrated parameters relevant for phytoplankton and zooplankton growth and turnover as described in Kriest et al. (2017) against observed global distributions of nutrients and oxygen.
Parameters to be calibrated for this new model were the light and nutrient affinities of phytoplankton, zooplankton quadratic 10 mortality, detritus remineralisation rate, particle stickiness and the exponent that relates particle sinking speed to particle size (see Table 2). After introduction of particle aggregation, the calibrated nutrient affinity of phytoplankton is now much higher, with a half-saturation constant for phosphate of K PHY = 0.11 mmol PO 4 m -3 instead of 0.5 mmol PO 4 m -3 in , very likely because the optimisation compensates for the higher export (and lower recycling) of phosphorus and nitrogen. Possibly for the same reason, detritus remineralisation rate in the optimised model is increased from 0.05 d -1 to 0. 25 15 d -1 . Light affinity of phytoplankton deviates less from the value in the model without particle aggregation, but the quadratic mortality of zooplankton is strongly reduced (1.6 (mmol P m -3 ) -1 instead of 4.55 (mmol P m -3 ) -1 ); the latter might be regarded as an attempt of the optimisation to reduce the export of organic matter from the euphotic zone. The two parameters that affect aggregation and particle sinking remained at moderate values of α = 0.42 and = 0.72, i.e. close to those applied in earlier model experiments with aggregation (e.g. Kriest, 2002). The residual cost function J RMSE of this pre-calibrated model 20 with aggregation was 0.472, i.e. lower than noAgg MIT2.8 (J RMSE = 0.529), but somewhat higher than achieved with a model version optimised against nutrient and oxygen concentrations , that resulted in a misfit of J RMSE = 0.439.
In the sensitivity experiment described below we will examine, whether this remaining misfit can be reduced even further, and evaluate the model sensitivity to changes in the parameters of this highly complex module.

Sensitivity experiments at coarse resolution (MIT2.8) 25
In the coarser model configuration of MOPS, MIT2.8, a first sensitivity study of 36 model simulations with different aggregation parameters was performed (see Table 1). We varied the values of four aggregation parameters, which control the rate of aggregation and the sinking behaviour of particles. The first parameter is the stickiness a, i.e. the probability that after collision two particles stick together, which was set to values of 0.2, 0.5 and 0.8, respectively. The second parameter is the maximum particle diameter for size dependent aggregation and sinking, D L , set to values of 1, 2 and 4 cm. A small value of speed of a primary particle with values of 0.7, 1.4 and 2.8 m d -1 . One effect of a small value of w 1 is that it reduces the loss of organic matter from surface layers, and thus has a direct effect on the recycling of nutrients at the surface. At the same time, it also affects the maximum possible sinking speed of the entire detritus pool. Finally, the exponent that relates particle sinking to diameter, , is set to values of either 0.62 and 1.17. A high represents dense particles, and a fast increase of particle sinking speed with size, a low value stands for more porous particles, which show only a weak relationship between 5 size and sinking speed (Kriest, 2002).

Sensitivity experiments at fine resolution (ECCO1.0)
The occurrence of aggregates, and their transport to the ocean interior, can furthermore depend on physical dynamics (e.g. Kiko et al., 2017). Therefore, in a second step, we repeated some of the experiments presented above in the finer resolution version ECCO1.0 to investigate possible improvements at higher resolution. In particular, we repeated all MIT2.8-10 simulations with = 0.62 in this finer resolution configuration. Additionally, we carried out three more simulations with = 1.17 but with the smallest D L = 1 cm to prevent particles from sinking through more than one box per time step (see Table   1). All simulations together lead to 30 model runs in the finer resolution configuration. To compare the ECCO1.0 simulations directly with results from MIT2.8, we re-gridded the result from ECCO1.0 simulations onto the coarser MIT2.8 grid. 15

Model Assessment and Diagnostics
Because observational data of particle flux are either limited with regard to space and time (e.g. Gehlen et al., 2006) or are combined with assumptions, that yield no clear patterns (Gehlen et al., 2006;Henson et al., 2012;McDonnell and Buesseler, 2010), this study restricts the model assessment to observations of nutrients and oxygen, in combination with the model fit to volume and location of oxygen minimum zones. 20

Root Mean Squared Error of Tracers
After a spin-up of 3,000 years into a seasonally cycling equilibrium state, the model results are evaluated in terms of annual means of oxygen, phosphate and nitrate. As in previous studies (e.g. Kriest et al., 2017) the misfit is calculated by the deviation between simulated results, m, and observed properties taken from the World Ocean Atlas (WOA), o, (Garcia et al., 2006). The deviations are weighted by volume of each grid box V i , expressed as the fraction of the total ocean volume V T . 25 The sum of the weighted deviations is normalised by the observed global mean concentration of each tracer: (1).
In this equation, j=1,2,3 describes the respective tracer (i.e. PO 4 , NO 3 and O 2 ). N is the total number of model grid boxes and o j is the global average observed concentration of each tracer . Thus, a low misfit value represents a good agreement between model and observations (J RMSE = 0 would be a perfect fit), which enables a prediction about the model 8 accuracy with regard to these tracers. The model runs with the lowest J RMSE in the coarse and the fine resolution are called RMSE MIT2.8* hereafter and RMSE ECCO1.0* , respectively.

Fit to oxygen minimum zones
To evaluate the extent and location of OMZs, we follow the approach of Cabré et al. (2015) by calculating the overlap between modelled and observed (Garcia et al., 2006; hereafter referred to as "WOA") OMZs. As several marine processes 5 are oxygen-dependent but have heterogeneous criteria for their minimum oxygen threshold, in this study, the OMZs are calculated for different oxygen threshold concentrations, C. Therefore, low-oxygen waters are characterised as O 2 < c, with c ranging from 0 to 100 mmol O 2 m -3 . To calculate the overlap between simulated and observed OMZs, we use the following equation (Sauerland et al., 2019): (2). 10 In this equation, ∩ is the volume of overlap of suboxic waters between model and observations, with regard to the defined oxygen threshold concentration c. This overlap is divided by the union (total volume of low-oxygen waters occupied in the model or in the observations) and results in a value between 0, equal to zero overlap between model and observations, and 1, which represents an optimal overlap. To adjust the scale to J RMSE , we calculated: In this equation, J OMZ varies between zero and one. Consequently, the scale of J OMZ is equivalent to the scale of J RMSE , which implies that a low misfit corresponds to a good agreement between model and observational data and vice versa. The model simulations with regard to lowest J OMZ are called OMZ MIT2.8* and OMZ ECCO1.0* hereafter. In calculating the overlap, we distinguish between the global ocean and the Pacific as well as the Atlantic Ocean. 20

Estimation of particle flux length scale b
To investigate, if, and how, the model reproduced observed maps of the particle flux length scale, b, that relates particle flux and depth via ∝ %& and derived from data by Marsay et al. (2015) and Guidi et al. (2015), we log-transformed F(z), the simulated, annual average flux of particulate organic matter as a function of depth and carried out a linear regression of these values. Highest b values correspond to short particle flux length scale, i.e. many small particles, and thus a low sinking 25 speed, shallow remineralisation and high oxygen consumption in shallow waters. For the reference models without aggregation these global maps should, in areas with shallow mixed layers, show spatially uniform values, as imposed by the model's prerequisites. Deviations from uniform values can either be ascribed to oxidant limitation of remineralisation (see above model description), or from physical processes such as mixing or upwelling, which can result in an additional vertical transport of particles.

9
The parameterisation of the aggregation model assumes a constant sinking speed for an upper size limit D L (see above), and therefore average particle sinking speed will remain constant below some depth. Also, the assumption of a particle size spectrum, size dependent sinking and constant remineralisation will result in particle flux profiles that do not fully agree with those predicted by a power law (see Kriest and Oschlies, 2008). Thus, because the aggregation model's prerequisites do not fully agree with a continuous increase of sinking speed with depth, we confine the regression of log-transformed particle flux 5 to a vertical range between 100-1,000 m, where the aggregation model still shows an increase of average sinking speed with depth (see also Kriest and Oschlies, 2008).

Global patterns of particle flux profiles
As could be expected, noAgg ECCO1.0 shows almost no spatial pattern of b, with values around the prescribed, nominal value 10 of b = 0.858 (global mean: 0.64; Fig. 1a; please note the different scaling in (a) and (d)) indicating long particle flux length scales and deep remineralisation. Regions with particularly low diagnosed b values (< 0.2) result either from decreased remineralisation in OMZs (e.g. eastern tropical Pacific OMZ) or are found in areas of deep mixing (in the model mainly high latitudes or western boundary currents), where vertical mixing increases the inferred particle flux length scales. However, for the best simulation with regard to the sum of J RMSE and J OMZ of the aggregation model (called ECCO1.0* hereafter) we find 15 highest b values, corresponding to short particle flux length scales, or shallow remineralisation, in the oligotrophic subtropical gyres. In contrast, b is smallest in the equatorial upwelling and in the shelf regions ( Fig. 1d and g). This pattern is in accordance with the observed spatial pattern derived by Marsay et al. (2015). In our model, this very deep flux penetration (b close to zero) in the equatorial upwelling can be explained with low oxygen concentrations, which reduce the remineralisation rate. In contrast, when deriving the particle flux length scale from a similar model but with oxygen-20 independent remineralisation (Kriest and Oschlies, 2013), we find a b close to the prescribed b value of 0.858 (Fig. S1).
In the subtropical and the equatorial region, the spatial variance (marked transparent red ; Fig 1g) of model-derived b values is quite high, which is caused by spatial variations in the physical environment, i.e. permanently stratified subtropical gyres and upwelling regions with low oxygen and reduced remineralisation. However, besides ECCO1.0* the four best model simulations with respect to the sum of J RMSE and J OMZ (simulation #14, #17, #28 and #29; Table 1) show essentially the same 25 pattern of b (Fig. S2), although these four simulations include quite different parameterisations (see Table 1).
Regions with high b values are characterised by a high spectral slope of the size distribution and therefore a high abundance of small particles, leading to slow sinking speeds (Fig. 7) and low export rates in ECCO1.0* (Fig. 1f). ECCO1.0* simulates highest export rates at high latitudes and in the upwelling region and lowest export rates in the subtropical gyres ( Fig. 1f and i). Although the spatial pattern of export rates is similar for both model simulations with and without aggregation, 30 ECCO1.0* shows a 1.6-fold higher global mean export rate (10.1 mmol P m -2 a -1 ) than noAgg ECCO1.0 (6.1 mmol P m -2 a -1 ). In ECCO1.0* export rates show a higher regional variability than in noAgg ECCO1.0 (Fig. 1c, 1f and 1i), which is due to blooms in the high latitudes during summer season accelerating the size-dependent aggregation and thus the export signal.
The oxygen concentration at a depth of 100 m shows the same global pattern in both simulations, with high oxygen concentrations at high latitudes, and decreasing concentrations towards the equator (Fig. 1b and 1e). However, the oxygen concentration at high latitudes is slightly higher in noAgg ECCO1.0 than in ECCO1.0* (Fig. 1h). Moreover, the global suboxic 5 volume (for a criterion c = 50 mmol m -3 ) in ECCO1.0* (7.3x10 16 m 3 ) is larger than in noAgg ECCO1.0 (3.7x10 16 m 3 ).
Comparing our model results with the dataset of Garcia et al. (2006), which yields a volume of 5.6x10 16 m 3 , we find an underestimation of the suboxic volume for noAgg ECCO1.0 by 34% and an overestimation for ECCO1.0* by 30%.

Representation of oxygen minimum zones
The finer resolution and data-assimilated circulation of ECCO1.0 in general improves the representation of OMZs in 10 comparison to MIT2.8 with regard to the overlap of OMZs for a criterion of 50 mmol m -3 (Fig. 2). Both simulations without explicit particle dynamics, namely noAgg MIT2.8 and noAgg ECCO1.0 , clearly underestimate the extent of the OMZ at a depth of 500 m and 1,000 m for an OMZ-criterion of 50 mmol m -3 in the Pacific basin (Fig. 2). The simulations including particle dynamics that are best with respect to the OMZ metric, OMZ MIT2.8* and OMZ ECCO1.0* , exhibit a larger OMZ area for both resolutions (Fig. 2). Despite the improved representation of OMZs, all models including the particle aggregation module still 15 tend to merge the OMZs of the Northern Hemisphere (NH) and the Southern Hemisphere (SH) at a depth of 500 m, which does not agree with the well separated northern and southern OMZ shown by the observations (Fig. 2 and Fig. S3). As reflected in a plot that shows the extent of OMZ in the northern and southern hemisphere, similar to Fig. 1a and 1b of Cabre et al. (2015), all models fail to represent the double structure of OMZ north and south of the equator. However, in our model the northern Pacific OMZ is fitted quite well ( Fig. 2 and Fig. S3). 20 Aggregation improves the representation of OMZs with respect to a criterion of c = 50 mmol m -3 compared to the simulations without aggregation for both resolutions in the NH, but not in the SH (Fig. 3). In noAgg ECCO1.0 the OMZ simulated in the NH is too small and too shallow (Fig. 3a). Even though OMZ ECCO1.0* tends to underestimate the suboxic area between ~700 m and 1,300 m, it shows a considerably higher overlap of model results and observations compared to noAgg ECCO1.0 (Fig. 3b). However, in the SH noAgg ECCO1.0 represents the OMZs better than OMZ ECCO1.0* , which tends to 25 overestimate the suboxic area in this hemisphere. In addition to differences caused by particle dynamics, circulation affects the performance in the two hemispheres: OMZ ECCO1.0* represents the highest overlap between ~100 and 500 m depth in the SH but this is surpassed by OMZ MIT2.8* between 500 and 900 m depth. In the NH, OMZ ECCO1.0* outcompetes OMZ MIT2.8* between 300 and 900 m depth as far as overlap is concerned (Fig. 3b).
However, the improvement of the representation of OMZs in the simulations with aggregation depends on the criterion for 30 OMZs. As could be expected, a higher oxygen threshold for the OMZ-criterion enhances the overlap between model simulations and observational data (Fig. 4). As for the fixed criterion of 50 mmol m -3 , globally and in the Pacific the better circulation and finer resolution of ECCO1.0 improves the overlap for varying OMZ-criterions in comparison to MIT2.8 ( Fig.  4a and c). While the OMZ ECCO1.0* simulation reaches globally a maximum overlap of 65.9% (for c = 100 mmol m -3 ), OMZ MIT2.8* culminates only in a maximum of 58.7% for the same criterion.
In the Pacific basin OMZ ECCO1.0* reaches an agreement with observations of 19.9% overlap for a criterion of 20 mmol m -3 (Fig. 4c). The overlap then increases strongly until the 100 mmol m -3 criterion (68.2%). It is noteworthy that globally and in the Pacific area noAgg ECCO1.0 outperforms all models for a criterion of 20 mmol m -3 , where it shows an agreement of almost 5 31%. The Atlantic basin shows an inverse trend (Fig. 4b): here, OMZ MIT2.8* represents the OMZ better than OMZ ECCO1.0* (26% and 12.2%, respectively, for a criterion of 70 mmol m -3 ). Further, in this region, the ECCO1.0 model that performs best with respect to RMSE (RMSE ECCO1.0* ) outperforms OMZ ECCO1.0* over the full range of criteria (Fig. 4b). Thus, there are large regional differences in the model's response to different circulations and particle dynamics. Because the dataset of observations used for comparison does not contain any concentrations below 30 mmol m -3 in the Atlantic, all models show 10 no overlap at all in this basin.
In summary, the improvement of model fit with regard to J OMZ depends not only on particle dynamics, but also on the definition of OMZs (i.e. the OMZ criterion c), the model resolution as well as the region considered (Fig. 2, Fig. 3, Fig. 4). Table 3 shows that in six cases out of nine (MIT2.8), a model that represents porous particles ( = 0.62) outperforms the 15 corresponding model with a sinking speed that describes rather dense, cell-like particles ( = 1.17). The same applies for the higher resolution (ECCO1.0), where in two cases out of three a porous parameterisation improves the fit with regard to J RMSE (see Table 1). Also, both J RMSE and J OMZ of the "dense" parameterisations are never among the best five models with respect to either metric (see Table 1). Thus, in the following we focus on model simulations with = 0.62.

Sensitivity of nutrient and oxygen distributions to aggregation parameters
Among the sensitivity experiments performed, the best model with respect to J RMSE (hereafter referred to as RMSE MIT2.8* ) is 20 characterised by an intermediate stickiness α of 0.5, the largest diameter for size-dependent aggregation and sinking, D L , of 4 cm and a minimum particle sinking speed w 1 of 2.8 m d -1, representing a rather fast organic matter transport to the ocean interior. However, many other models with medium stickiness perform about equally well (Fig. 5, upper mid panel). Models with lower stickiness perform best with slow minimum sinking speed w 1 and a large maximum size D L =4 cm for sizedependent sinking and aggregation (Fig. 5, upper left panel). In contrast, a large stickiness (which facilitates the formation of 25 aggregates in surface layers) requires either small w 1 or D L , which reduces the export of particles out of the euphotic zone, and into the ocean interior.
Oxygen concentrations contribute most to the global J RMSE . The influence of oxygen on global tracer misfit is dominated by the deep concentrations (Fig. S4), and thus to a large extent by the large-scale circulation. The OMZs, because of their small regional extent, contribute less to the global misfit . This is confirmed by Fig. S4  30 (d, e, f), showing that, in the eastern tropical Pacific region deep (>300 m) mesopelagic and deep oxygen concentrations scatter strongly among the different models (Fig. S4 a), despite their good global match in shallow waters. Likewise, although global mean profiles of nutrients are quite similar among the different circulations, and agree quite well with observations, their concentrations scatter strongly in the eastern tropical Pacific. Most of the simulations tend to underestimate the oxygen and nitrate concentration in this region (Fig. S4 a and c). Too low oxygen concentrations lead to too high denitrification and thus widespread nitrate depletion in the eastern tropical Pacific region, which explains the simultaneous underestimate of oxidants in this region.
To sum up, a moderate stickiness enhances the chance of a good model fit to nutrients and oxygen (J RMSE ), but there is no 5 unique trend for the parameters or combination of parameters, with the exception of the exponent that relates particle sinking speed to its size: here, we find an advantage of a parameterisation characteristic for porous marine aggregates. In the optimal scenario, the misfit is less than that of a model without aggregates, when this is simulated with fixed reference parameters (noAgg MIT2.8 ). Because of the small spatial extent of OMZs, the model fit to nutrient and oxygen concentrations is mainly caused by the large-scale tracer distribution, even if some models show a considerable mismatch to these tracers in OMZs. 10 The pattern for J RMSE does not change very much when applying a different, higher resolved and data assimilated circulation (see Table 1 and Fig. 6). Now, the optimal model (RMSE ECCO1.0* ) is improved with respect to J RMSE by about 13%, but many other, almost equally good solutions, can be found with moderate to high stickiness. Introducing aggregates in this coupled model system does not improve the model fit to nutrient and tracer concentrations, as evident from the comparison of RMSE ECCO1.0* (J RMSE = 0.431) against a model without aggregate dynamics (J RMSE = 0.426; Table 1). The lack of 15 improvement can likely be explained by the fact that the biogeochemical parameters of MOPS with particle dynamics were adjusted in the circulation of MIT2.8, and thus are not optimal for the model when simulated in the physical dynamics of ECCO1.0.
The sensitivity to the metric for OMZs differs from the one to the metric for nutrients and oxygen. Now, for the fit to oxygen minimum zones (J OMZ ), a large stickiness, α, in combination with D L of 2 cm and slow to moderate minimum sinking speed 20 w 1 are of advantage ( Fig. 5 and Fig. 6). Thus, a high rate of aggregation, and a maximum sinking speed of about 50-100 m d -1 improves the model with respect to OMZs. This is also evident from comparison of the optimal models (OMZ MIT2.8* and OMZ ECCO1.0* ) to models without aggregate dynamics (noAgg MIT2.8 and noAgg ECCO1.0 ), shown in Fig. 3 and Fig. 4 and subsection 3.2. Nevertheless, even the models that perform best with respect to J OMZ underestimate mesopelagic oxygen when averaged over the eastern tropical Pacific (Fig. S4 a). 25 The sensitivity patterns with regard to J OMZ among both configurations MIT2.8 and ECCO1.0 diverge considerably from each other, which is in contrast to the patterns for J RMSE noted above (compare Fig. 5 with Fig. 6). Thus, model performance with respect to J OMZ seems to depend much more on circulation and physical details than the large-scale dynamics reflected in J RMSE .

Discussion
aggregates, which are composed of phytoplankton and detritus, the parameterisation, which is based on dense particles (dSAM, Kriest 2002) and a biogeochemical model, which is different. We found high values for the spectral slope of the size distribution (i.e. high abundance of small particles) and thus a low particle sinking speed in the subtropical gyres (Fig. 7), which corresponds with the findings by Oschlies and Kähler (2004) and Dutay et al. (2015). This, in turn, leads to highest b values in the oligotrophic subtropical gyres and lowest ones in the high latitudes and the upwelling region, and agrees with 5 the pattern as shown in Marsay et al. (2015). These findings imply that such a b pattern cannot only result from temperature dependent remineralisation -as suggested by Marsay et al. (2015) -but also from particle dynamics and temperatureindependent remineralisation. However, if temperature-dependent remineralisation, as suggested by Marsay et al. (2015) or Iversen and Ploug (2013), was also included in our model, this would likely enhance horizontal variations in the particle flux profile, with even deeper flux penetration in the cold waters of the high latitudes and upwelling areas. Beside particle 10 dynamics, the low b values in upwelling regions found in our study (Fig. 1d), are also caused by the suboxic conditions, which suppress remineralisation in subsurface waters. Such a tight link between suboxia and deep flux penetration is supported by the observations reported by Devol and Hartnett (2001) and Van Mooy et al. (2002). Therefore, two different processes -particle aggregation and/or temperature-dependent remineralisation -suggest low b values and deep flux penetration in the very productive areas of high latitudes. A third process, which consists in oxygen-dependent 15 remineralisation, is superimposed on these in OMZs, causing the steepest particle profiles in these areas.
However, it should be noted that although the maximum sinking speed of our best simulations (101 m d -1 (#17) and 51 m d -1 (#26), see Table 1) agrees with observations (Alldredge and Gotschalk, 1988;Nowald et al., 2009, Jouandet et al., 2011, the range of b values in our model is almost twice as large as suggested by most empirical studies (Berelson, 2001;Buesseler et al., 2007;Martin et al., 1987;Van Mooy et al., 2002). However, as there is no common depth range to determine the particle 20 flux length scale b, the depth range spreads over a wide range in various studies and thus impedes the comparability (Marsay et al. 2015), which might explain some divergence between observations and model results. In particular, our model simulates a too large fraction of small particles and therefore a too steep particle size spectrum in the subtropical gyres, which causes too high b values in these areas. Other processes that modify the size spectrum, like grazing by zooplankton, and the subsequent egestion of large fecal pellets, might also play a role in these regions. Additionally, the model tends to 25 underestimate the number of large particles (size range 0.14 to 16.88 mm) in the surface of the tropical Atlantic Ocean (23°W), compared to observations (Kiko et al., 2017;Fig. S6). On the other hand, a first, direct comparison to the UVP 5 dataset (Kiko et al., 2017, their Fig. 1) exhibits a correct magnitude regarding the number of particles within this size range (0.14 to 16.88 mm) in our model ( Figure S5) along the 151°W section. One possible explanation for the mismatch at 23°W could consist in a not sufficiently resolved equatorial current system, which also will be discussed below. Also, additional 30 biological processes such as the downward transport of organic matter trough vertically migrating zooplankton (Kiko et al., 2017), or particle breakup of aged, fragile particles at depth (e.g. Biddanda et al., 1988) could improve the model. However, introducing this additional complexity is beyond the scope of this paper. In future studies, consideration of these processes, in conjunction with a comprehensive model calibration against observed particle abundances and size spectra (e.g.
contributions of individual processes such as aggregation, vertical migration, temperature-dependent remineralisation and to validate simulated particle dynamics.
However, model calibration against observed particle dynamics has to account for characteristics and limitations of observations. For example, the size spectrum assumed in our model is of infinite upper size and also contains particles with a 5 diameter larger than, e.g. 4 cm (the upper limit for size-dependency of aggregation and sinking). While these particles exist (e.g. Bochdansky and Herndl, 1992), they are very rare (in the model, and likely also in the observations), and might not be observed with standard methods, which usually rely on a sample size of few litres. The rare occurrence of large particles, and the limited sample size has, for example, consequences for estimated size spectra parameters (Blanco et al., 1994). Thus, any model calibration against observations of particle abundance and size has to account for a proper match between simulated 10 and observed quantities.
As we used on the one hand two different model grid resolutions and on the other hand varied model parameterisations with regard to particle aggregation, changes in the location and extension of OMZs and the distribution of tracers within each resolution are exclusively driven by the aggregation parameters. A good parameterisation of particle aggregation parameters can therefore have a major influence on the representation of OMZs. Furthermore, a higher model resolution improves the 15 depiction of equatorial currents and therefore the oxygen transport (Cabré et al., 2015;Duteil et al., 2014), which, in turn, results in an improved representation of OMZs in the finer resolution configuration, ECCO1.0, compared to the coarser resolution, MIT2.8. However, as physical processes at smaller scales affect the simulated shallow to mesopelagic oxygen and nutrient concentrations for the eastern tropical Pacific (Getzlaff and Dietze, 2013), the finer (1°x1°) resolution of ECCO1.0 is not sufficient to resolve the details of the equatorial current system (Duteil et al., 2014). This can explain the 20 still high residual misfit of these simulations, and the missing double structure of OMZs in the Eastern Tropical Pacific. We therefore suggest, that the difference in improving the representation of OMZs between northern and southern hemisphere is more affected by physics than by biology.
Furthermore, results of our sensitivity study confirm that dense particles do not constitute a realistic representation of particles, as indicated by Karakaş et al. (2009) and Kriest (2002). Porous particles seem to constitute a more appropriate 25 parameterisation for good model fit with regard J RMSE and J OMZ (Table 1). Although the observed stickiness ranges between almost zero and one (e.g. Alldredge and McGillivary, 1991;Kiørboe et al., 1990), in our study a moderate stickiness, a, between 0.5 and 0.8 leads the model towards a good fit to observed nutrients, oxygen and OMZs.
In summary, our study supports the results of Schwinger et al. (2016), who found an improved representation of nutrient distribution and OMZs when switching from constant particle sinking to either a power law or particle dynamics, similar to 30 those presented here. However, the difference between the two latter schemes in that study were only small. A more extensive search of the parameter space within a given circulation may further improve the model. Additionally, we optimised noAgg MIT2.8 against the same misfit function as MOPS oD of Kriest et al. 2017 and found that even though including an aggregation module improves our model, utilising an appropriate parameter optimisation would further enhance our model fit. Thus, without a comprehensive calibration of biogeochemical and aggregation parameters there only seems to be a slight advantage when using this more complex model of particle dynamics.
Finally, we found a steep particle size spectrum in the subtropical oligotrophic region (Fig. 1d), which does not agree with observational data. Potentially, there are processes taking place, which are not considered in our model i.e. particle repackaging and active transport by zooplankton (vertical migration) (Kiko et al. 2017) based on a modified food web. Thus, 5 particle aggregation alone so far seems not to be sufficient for a correct representation of the particle size spectrum. Najjar et al. (2007) applied different model circulations to the same biogeochemical model, and found that that physical processes are an important factor for modelling marine biogeochemistry. Our study furthermore showed that also biogeochemical parameterisations -in particular, those related to particle flux -can have an important impact on the 10 representation of dissolved inorganic tracers, in line with earlier studies (e.g. Kriest et al., 2012;Primeau, 2006, 2008). These earlier studies applied and varied a globally uniform particle flux length scale, whereas it has been suggested that this parameter should vary in space and time (e.g. Guidi et al., 2015;Marsay et al., 2015). The sensitivity study presented here constitutes a first approach to systematically estimate the impact of marine particle aggregation -and thus a spatially and temporally variable flux length scale -on the location and extent of OMZs as well as the representation of 15 phosphate, nitrate and oxygen under steady-state conditions in a global three-dimensional biogeochemical ocean model.

Conclusion and Outlook
We have shown that the assumptions inherent in the model confirm the general pattern of the spatial map of b values proposed by Marsay et al. (2015) (Fig. 1a and d). This, in turn, shows that the pattern of Martin's b cannot only be depicted by a POC flux dependent on temperature but also by simulating explicit particle dynamics.
We furthermore found that even though there are still a lot of gaps in understanding several processes e.g. the variation of 20 export rates, particle stickiness and particle flux profile over space and time, as well as the link between particle diameter and sinking speed, the comparisons against observational data show a trend towards a model improvement by integrating particle dynamics (Table 1). While the parameterisation of aggregation leads the model towards an improved fit to OMZs for both model resolutions, this increase in model fit with regard to phosphate, nitrate and oxygen is only detectable in the coarse resolution MIT2.8, but not in the finer resolution and data-assimilated circulation of ECCO1.0. Moreover, model 25 simulations show that besides effects of grid resolution, the model fit with regard to J RMSE and J OMZ is mainly driven by the particles' porosity. Our results indicate that a best fit to both, tracers as well as OMZs (50 mmol O 2 m -3 criterion), is achieved by parameterising porous particles in combination with an intermediate to large maximum particle diameter for size dependent aggregation and sinking, a moderate to high stickiness ranging between 0.5 and 0.8 and an intermediate to high initial sinking speed ranging between 1.4 and 2.8 m d -1 (Fig. 5). The strong sensitivity of the model fit to aggregation 30 parameters may point towards the importance of a spatially and temporally varying flux length scale; however, they also show, that the dynamics of the model depend strongly on the assumptions we make with respect to particle properties and processes.
Finally, we have shown that uncertainties in the parameterisation of particle aggregation remain, leading to the inference that dissolved inorganic tracers offer only insufficient observational constraints for global particle parameterisation. Therefore, for an accurate representation it will be necessary to calibrate the model not only against observed phosphate, nitrate, oxygen 5 distributions and volume and location of OMZs (Sauerland et al., 2019) but also against number and size of particles, using comprehensive datasets of observations (as in Guidi et al., 2015).

Code and data availability
The source code of MOPS including the aggregation module coupled to TMM as well as the model output are available at: https://data.geomar.de/thredds/catalog/open_access/niemeyer_et_al_2019_bg/catalog.html. 10

Author contributions
D. Niemeyer, I. Kriest and A. Oschlies conceived the study. D. Niemeyer performed and analysed the simulations. All authors discussed and wrote the manuscript. Table 1: Model runs of sensitivity study, their parameter combinations and the calculated misfit of tracers (J RMSE ) and OMZs (J OMZ ) for MIT2.8 and the ECCO1.0 configurations. The 25% best simulations with regard to J RMSE and J OMZ are highlighted in yellow and the worst 25% in red (relative to RMSE ECCO1.0* and OMZ ECCO1.0* ). The simulations in between are coloured in two