2.1 Oceanic transport

In this study, we used the “transport matrix method” (TMM) (Khatiwala et al., 2005; Khatiwala, 2007, 2018), 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-resolution version, based on a data-assimilated circulation (ECCO1.0) (Stammer et al., 2004). The MIT2.8 forcing and transport represent a resolution of $\mathrm{2.8}{}^{\circ }\times \mathrm{2.8}{}^{\circ }$ and 15 depth layers with a thickness ranging between 50 and 690 m. ECCO1.0 TMs and forcing are based on a resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ and 23 depth layers, with a thickness ranging between 10 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 d for physical transport and a time step length of 1∕16 d for biogeochemical interactions in the coarse resolution, MIT2.8. Because some parameter configurations allow a very large 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 d 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 d (see Table 1) but with the exception of three experiments, where we used a length of 1∕160 d (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 3000 years until tracers approached steady state. The last year is used for analysis as well as misfit calculations.

Table 1Model 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${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ and OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$). The simulations in between are coloured in two orange gradations (bright orange is medium good and dark orange is medium bad). The best simulation of each resolution with regard to J_RMSE and J_OMZ is bold. OMZ is defined as 50 mmol m⁻³. Parameter η denotes the exponent for size-dependent sinking, α the stickiness, w₁ the minimum sinking speed, D_L the maximum diameter for size-dependent sinking and aggregation, and w_max the maximum sinking velocity in the spectral computations.

Download Print Version

2.2 The biogeochemical model

2.3 Model simulations and experiments

2.3.2 Adjustment of biogeochemical model parameters

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 Kriest (2017), 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 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₄ m⁻³ instead of 0.5 mmol PO₄ m⁻³ in Kriest et al. (2017), 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 to 0.25 d⁻¹. 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${}^{-\mathrm{3}}{)}^{-\mathrm{1}}$ instead of 4.55 (mmol P m${}^{-\mathrm{3}}{)}^{-\mathrm{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 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 (Kriest et al., 2017), which 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.

Table 2Model adjustment of biogeochemistry with aggregates compared to Kriest et al. (2017) and new parameters in this study.

Download Print Version | Download XLSX

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

3.1 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 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 the highest b values, corresponding to short particle flux length scales, or shallow remineralisation, in the oligotrophic subtropical gyres. In contrast, b is the 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 0) 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-independent remineralisation (Kriest and Oschlies, 2013), we find a b close to the prescribed b value of 0.858 (Fig. S1 in the Supplement).

Figure 1Global maps of b (a, d), O₂ at 100 m (mmol m⁻², b, e) and export at 100 m (mmol P m⁻² a⁻¹, c, f) for noAgg^ECCO1.0 (a, b, c) and for the best aggregation model with regard to the sum of J_RMSE and J_OMZ (simulation no. 26; d, e, f). The black line indicates the OMZ for a criterion of 50 mmol m⁻³. Lower panels: Global mean (dotted line) and standard deviation (transparent shaded) of b (g), O₂ (h) and export (i) of noAgg^ECCO1.0 (black) and the best aggregation model with regard to the sum of J_RMSE and J_OMZ (simulation no. 26; red). Please note the different scaling for b values (a, d).

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 nos. 14, 17, 28 and 29; Table 1) show essentially the same 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 the highest export rates at high latitudes and in the upwelling region and the 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, ECCO1.0^* shows a 1.6-fold higher global mean export rate (10.1 mmol P m⁻² a⁻¹) than noAgg^ECCO1.0 (6.1 mmol P m⁻² a⁻¹). In ECCO1.0^* export rates show a higher regional variability than in noAgg^ECCO1.0 (Fig. 1c, f and i), 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 e). However, the oxygen concentration at high latitudes is slightly higher in noAgg^ECCO1.0 than in ECCO1.0^* (Fig. 1h). Moreover, the global suboxic volume (for a criterion c=50 mmol m⁻³) in ECCO1.0^* (7.3×10¹⁶ m³) is larger than in noAgg^ECCO1.0 (3.7×10¹⁶ m³). Comparing our model results with the dataset of Garcia et al. (2006), which yields a volume of 5.6×10¹⁶ m³, we find an underestimation of the suboxic volume for noAgg^ECCO1.0 by 34 % and an overestimation for ECCO1.0^* by 30 %.

3.2 Representation of oxygen minimum zones

The finer-resolution and data-assimilated circulation of ECCO1.0 in general improves the representation of OMZs in comparison to MIT2.8 with regard to the overlap of OMZs for a criterion of 50 mmol m⁻³ (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 and 1000 m for an OMZ criterion of 50 mmol m⁻³ in the Pacific basin (Fig. 2). The simulations including particle dynamics that are the best with respect to the OMZ metric, OMZ${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$ and OMZ${}^{\mathrm{ECCO}{\mathrm{1.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 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 (Figs. 2 and S3 in the Supplement). As reflected in a plot that shows the extent of OMZ in the NH and SH, similar to Fig. 1a and 1b of Cabré 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 (Figs. 2 and S3).

Figure 2Comparison of Pacific Ocean OMZ (O₂ ≤ 50 mmol m⁻³) between model simulations and observations. Panels (a) and (b) show the OMZ at a depth of 500 and 1000 m for the coarse resolution, MIT2.8, and panels (c) and (d) for the fine resolution, ECCO1.0.

Aggregation improves the representation of OMZs with respect to a criterion of c=50 mmol m⁻³ 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${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ tends to underestimate the suboxic area between ∼700 and 1300 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${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$, which tends to overestimate the suboxic area in this hemisphere. In addition to differences caused by particle dynamics, circulation affects the performance in the two hemispheres: OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ represents the highest overlap between ∼100 and 500 m depth in the SH, but this is surpassed by OMZ${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$ between 500 and 900 m depth. In the NH, OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ outcompetes OMZ${}^{\mathrm{MIT}{\mathrm{2.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 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⁻³, globally and in the Pacific the better circulation and finer resolution of ECCO1.0 improves the overlap for varying OMZ criteria in comparison to MIT2.8 (Fig. 4a and c). While the OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ simulation reaches globally a maximum overlap of 65.9 % (for c=100 mmol m⁻³), OMZ${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$ culminates only in a maximum of 58.7 % for the same criterion.

In the Pacific basin OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ reaches an agreement with observations of 19.9 % overlap for a criterion of 20 mmol m⁻³ (Fig. 4c). The overlap then increases strongly until the 100 mmol m⁻³ 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⁻³, where it shows an agreement of almost 31 %. The Atlantic basin shows an inverse trend (Fig. 4b): here, OMZ${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$ represents the OMZ better than OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ (26 % and 12.2 %, respectively, for a criterion of 70 mmol m⁻³). Further, in this region, the ECCO1.0 model that performs best with respect to RMSE (RMSE${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$) outperforms OMZ${}^{\mathrm{ECCO}{\mathrm{1.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⁻³ in the Atlantic, all models show 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 (Figs. 2, 3, 4).

Table 3Number of simulations with different parameters for D_L, α and w₁ for the porous (η=0.62) and dense (η=1.17) particles which outperform the corresponding other size. The numbers are given with respect to two different criteria, J_RMSE and J_OMZ.

Download Print Version | Download XLSX

3.3 Sensitivity of nutrient and oxygen distributions to aggregation parameters

Table 3 shows that in six cases out of nine (MIT2.8), a model that represents porous particles (η=0.62) outperforms the 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.

Among the sensitivity experiments performed, the best model with respect to J_RMSE (hereafter referred to as RMSE${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$) is 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₁ of 2.8 m d⁻¹, representing a rather fast organic matter transport to the ocean interior. However, many other models with medium stickiness perform about equally well (Fig. 5b). Models with lower stickiness perform best with slow minimum sinking speed w₁ and a large maximum size D_L=4 cm for size-dependent sinking and aggregation (Fig. 5a). In contrast, a large stickiness (which facilitates the formation of aggregates in surface layers) requires either small w₁ 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 (Kriest et al., 2017). 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 (Kriest et al., 2017). This is confirmed by Fig. S4d, e and f, showing that, in the eastern tropical Pacific region, deep (> 300 m) mesopelagic and deep oxygen concentrations scatter strongly among the different models (Fig. S4a), 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. S4a and c). Oxygen concentrations that are too low lead to denitrification that is too high 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 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.

The pattern for J_RMSE does not change very much when applying a different, more highly resolved and data-assimilated circulation (see Table 1 and Fig. 6). Now, the optimal model (RMSE${}^{\mathrm{ECCO}{\mathrm{1.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${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$ (J_RMSE=0.431) against a model without aggregate dynamics (J_RMSE=0.426; Table 1). The lack of 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 they are not optimal for the model when simulated in the physical dynamics of ECCO1.0.

Figure 3Area of OMZ (a) and overlap of OMZs between model and observations following Cabre et al. (2015) (b). In both panels, the left-hand side shows the Southern Hemisphere (0–40^∘ S); the right shows the Northern Hemisphere (0–40^∘ N), plotted against the logarithmic depth. OMZs are defined as regions with O₂ < 50 mmol m⁻³.

Download

Figure 4Overlap between modelled and observed OMZs (Eq. 2) for varying criteria c, ranging from 0<c<100 mmol m⁻³) on a global scale (a), for the Atlantic Ocean (b) and for the Pacific Ocean (c).

Download

The sensitivity to the metric for OMZs differs from the sensitivity 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 w₁, is of advantage (Figs. 5 and 6). Thus, a high rate of aggregation, and a maximum sinking speed of about 50–100 m d⁻¹, improves the model with respect to OMZs. This is also evident from comparison of the optimal models (OMZ${}^{\mathrm{MIT}{\mathrm{2.8}}^{*}}$ and OMZ${}^{\mathrm{ECCO}{\mathrm{1.0}}^{*}}$) to models without aggregate dynamics (noAgg^MIT2.8 and noAgg^ECCO1.0), shown in Figs. 3 and 4 and Sect. 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. S4a).

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.

Figure 5Sensitivity of J_RMSE (Eq. 1); (a, b, c) and J_OMZ (Eq. 3); (d, e, f) to minimum sinking speed w₁ and maximum size D_L for the coarse resolution, MIT2.8, for three different values of stickiness (a–f), and η=0.62 (“porous” particles). The colour bar shows J_RMSE and J_OMZ (blue – good fit, red – bad fit), normalised by its minimum value across all model experiments. Black arrows indicate an improvement of J_RMSE or J_OMZ with increasing parameter values, while white arrows show an improvement with decreasing values.

Download

Figure 6As Fig. 5 but for simulations with ECCO1.0.

Download