A new mathematical model to explore microbial processes and their constraints in phytoplankton colonies and sinking marine aggregates

See allHide authors and affiliations

Science Advances  31 Oct 2018:
Vol. 4, no. 10, eaat1991
DOI: 10.1126/sciadv.aat1991


N2-fixing colonies of cyanobacteria and aggregates of phytoplankton and detritus sinking hundreds of meters per day are instrumental for the ocean’s sequestration of CO2 from the atmosphere. Understanding of small-scale microbial processes associated with phytoplankton colonies and aggregates is therefore crucial for understanding large-scale biogeochemical processes in the ocean. Phytoplankton colonies and sinking aggregates are characterized by steep concentration gradients of gases and nutrients in their interior. Here, we present a mechanistic mathematical model designed to perform modeling of small-scale fluxes and evaluate the physical, chemical, and biological constraints of processes that co-occur in phytoplankton colonies and sinking porous aggregates. The model accurately reproduced empirical measurements of O2 concentrations and fluxes measured in sinking aggregates. Common theoretical assumptions of either constant concentration or constant flux over the entire surface did not apply to sinking aggregates. Consequently, previous theoretical models overestimate O2 flux in these aggregates by as high as 15-fold.


The rising CO2 concentration in the Earth’s atmosphere is partly counterbalanced by CO2 sequestration in phytoplankton and N2-fixing cyanobacteria. One of the major topics in aquatic microbial ecology deals with the fate of fixed CO2 in marine aggregates composed of phytoplankton and detritus and the pathways by which they are transported downward through the water column, decomposed and mineralized by microbes, and grazed by zooplankton, that is, the ocean’s biological carbon pump (1, 2). The settling of marine aggregates in the water column has substantial effects on the ecology and the biogeochemistry of carbon and nitrogen of the oceans (3). Because of their relatively high sinking speeds, aggregates might be thought of as vehicles for vertical flux of organic matter in aquatic environments. However, marine aggregates also constitute hot spots of microbial activity and are sites of rapid and efficient turnover of particulate organic carbon in the ocean (1). Solubilization and remineralization of organic substrates often imply that concentrations of oxygen, nutrient, and dissolved organic carbon within aggregates are significantly different from those in the bulk water (46), and the measured influx of oxygen and efflux of dissolved organic matter (DOM) such as combined amino acids from aggregates suggest efficient mass transfer across the aggregate-water interface (79). The quantification of exchange of solute fluxes between the aggregate and the surrounding water is still an important topic of ongoing research. Estimates of total O2 uptake by individual aggregates are usually made either by laboratory experiments or by mathematical models based on mass transfer theories. Experimental determination of O2 uptake is often based on calculation of O2 flux using O2 concentration profiles measured with microsensors. In this method, the microsensor approaches an immersed aggregate and measures the concentration profile, usually along the downstream axis. The concentration profile is then used to calculate the local flux. Hereby, the aggregate is stabilized in the suspended phase via an upward-directed flow velocity that opposes and balances the sinking velocity of the aggregate (8, 9). By considering a uniform concentration gradient around the aggregate, the concentration gradient measured at downstream pole (perpendicular to the aggregate’s surface) is then multiplied by diffusion coefficient and the aggregate’s surface area to obtain the total uptake (2, 8, 10). In an attempt to quantify the O2 uptake of sinking marine aggregates, Kiørboe et al. (11) considered sinking marine aggregates as impermeable to flow, and the fluxes at the aggregate-water interface were modeled, assuming constant concentration (Dirichlet boundary condition) or constant flux (Neumann boundary condition) at the aggregate’s surface to solve the advection-diffusion equation. The Dirichlet or Neumann boundary condition was used, assuming that the solute transfer was either transport limited or reaction limited. It should be noted that the model by Kiørboe et al. (11) had no reaction term; only mass transfer at the aggregate-water interface was considered. However, it is well known that most of marine aggregates are highly porous, with porosities reaching up to 99% (12). Thus, considering aggregates as solid spheres in models may not adequately represent in situ aggregates. In addition, when model aggregates are assumed to be solid (11), it is implied that bacterial activities can only occur at the aggregate’s surfaces, contradicting the observations. Moreover, when addressing specific biological processes such as O2 respiration by living organisms, including bacteria, phytoplankton, and protozoans, within the aggregates, the advection-diffusion equation has to be extended by inclusion of an appropriate reaction term, for example, production or consumption of gases, solutes, and/or nutrients inside aggregates. A model for mass transfer within and around sinking porous aggregates was already developed by Liu et al. (13) but without a reaction term. To overcome the aforementioned limitations, we designed a model for porous, impermeable aggregates that is capable of accounting for biological activities within the aggregates by adding a reaction term. The model can also be used to address physical, chemical, and biological constraints and rates of various processes that may (co-)occur in the microenvironment of sinking aggregates. As an example, we used this model to simulate the entire flow and O2 field within and around sinking marine aggregates and accurately quantify total O2 uptake by microbes associated with the aggregates.


Depending on the size and excess density of marine aggregates, their sinking velocities fall in a range from 10 m/day to several hundred meters per day (2, 12). For an aggregate of radius R sinking with the velocity U0, the corresponding Reynolds number reads Re = U0R/υ, with υ being the kinematic viscosity of the ambient water. Empirical measurements reveal that Re for marine aggregates often is in the range between 0.1 and 10 (12). In the present study, stagnant aggregates subjected to a constant upward velocity U in the far field (UU0) were considered. The geometrical configuration is sketched in Fig. 1A. To calculate the flow field around sinking aggregates, we solved the corresponding Navier-Stokes equations (see Eq. 5 and Materials and Methods for flow field). Although the flow fields for a wider Re range were calculated, we here display the streamlines for Re = 0.1 and Re = 10 only (Fig. 1, B and C). In all simulations, the aggregate’s radius and fluid viscosity were kept constant, and therefore, a change in Reynolds number is due to change of the velocity only. The flow is almost symmetrical upstream and downstream for the case of Re = 0.1 (Fig. 1B). As Re increases, the upstream-downstream symmetry breaks, and the streamlines are displaced by the aggregate for a larger distance behind it (downstream) than in front of it (upstream) (Fig. 1C). For the Re range considered here, this effect became more pronounced with increasing Re. For each Reynolds number, a region of slower fluid velocity (with respect to U) formed in the downstream, which extended with increasing Re. This can be better illustrated when plotting the distances between the downstream pole and those points above the aggregate at which the downstream velocity reaches 10% of the corresponding U, represented by ZU10, as a function of Re (Fig. 1D).

Fig. 1 Definitions of model parameters, simulated velocity fields around sinking aggregates, and sensitivity analyses.

(A) Sketch of the modeled sinking aggregates. (B and C) Stream lines around sinking aggregates at Re = 0.1 and Re = 10, respectively. The color bars show the magnitude of fluid velocity normalized by the far-field velocity U for each case. (D) The distance from the downstream pole at which the velocity reaches 10% of U along the z axis. (E to G) The O2 concentration profiles along the z axis for different Da (Re = 1.4, Km = 1, ε = 0.95), different Km (Re = 1.4, Da = 1, ε = 0.95), and different ε (Re = 1.4, Da = 1, Km = 1), respectively.

Once the flow field is calculated, one can obtain the concentration field of any solute inside and around the aggregate by solving the concentration equation (see Eq. 6 and Materials and Methods for concentration field). The oxygen concentration field within and around porous aggregates depends on the ambient concentration, aggregate porosity, Reynolds number, and the reaction term. For the computation of concentration fields at a given Reynolds number, we set the O2 diffusion coefficient, ambient concentration, porosity, half-saturation constant, and Damköhler number to D = 1.70 × 10−5 cm2/s, C = 300 μM, ε = 0.95, Km = 1 μM, and Da = 1, respectively (see Materials and Methods for the concentration field). Values for ambient concentration and porosity are typical and have been reported in the literature frequently. The choice of Km = 1 μM was based on estimates in previous studies (14). The rationale for taking Da = 1 was to best reproduce the measurements by Ploug (10), in which the O2 concentration at the aggregate’s center reached almost 82% of the ambient concentration at Reynolds numbers around 1 (see Fig. 1E and Materials and Methods for model validation). The O2 content inside the aggregate must decrease with an increase of microbial respiration. This fact has been shown in Fig. 1E. Note that the maximum reaction rate Vm is proportional to Da. For the case of Da = 10, the microbial O2 respiration rate is much larger than the diffusive rate, and thus, the O2 uptake by the aggregate is diffusion limited (Fig.1E). Furthermore, it is interesting to examine the effect of Km and ε variation on the O2 concentration profile when other parameters remain unchanged. Although the O2 concentration increases with the increase of Km, the concentration variation is not significant as long as Km is sufficiently lower than C (Fig. 1F). Since the interstitial water inside the aggregate is reduced with the decrease of porosity, the O2 concentration inside the aggregate decreases too (Fig. 1G). Figure 2 illustrates the normalized O2 concentration field, C/C, for different Reynolds numbers. For the case of Re = 0, that is, when advection was absent, the concentration field was spherically symmetric, as expected. Here, the normalized ambient O2 concentration of C/C = 1.0 falls within a distance of approximately 3R to 0.8 at the aggregate’s surface and to 0.65 at the aggregate’s center. With increasing Re, the concentration field around the aggregate elongated in the direction of flow, producing downstream plumes in the rear of aggregate. At higher Reynolds numbers, the plume behind the aggregate became longer and thinner. The response of increasing Reynolds number to O2 concentration profiles along the x (horizontal) and z (vertical) axes is shown in Fig. 3. Those computations revealed that the O2 concentration increased monotonically along the x axis and asymptotically approached the ambient concentration outside the aggregate. This trend was valid for the entire Re range examined here (Fig. 3A). In addition, O2 concentration increased at both the center and the surface of the aggregate when Re increased (Fig. 3, A and B). As shown, the distances at which the O2 concentrations closely approached the ambient concentration decreased with increasing flow. The O2 concentration profiles, however, become more compressed with increasing Re (for example, compare Re = 5 and Re = 10).

Fig. 2 Effect of sinking velocity on the concentration field.

O2 concentration field around and within sinking porous aggregates at Re = 0 (A), Re = 0.01 (B), Re = 0.05 (C), Re = 0.1 (D), Re = 1 (E), Re = 5 (F), and Re = 10 (G).

Fig. 3 Effect of sinking velocity on the concentration profiles.

(A) O2 concentration profiles along equatorial axis for different Re. The vertical dashed line shows the aggregate’s interface along the x axis. Because of the rotational symmetry around the z axis, the profiles have only been shown for x ≥ 0. (B) Zoom in on the interface along the equatorial axis. (C) O2 concentration profiles along the z axis for different Re. The left and right vertical dashed lines show the upstream and downstream interfaces, respectively. Zoom in on the upstream (D) and downstream (E) interfaces along the z axis.

Unlike O2 concentration along the x axis, concentration along the z axis was asymmetric (Fig. 3C). While the difference between the upstream O2 concentration profiles (Fig. 3D) and those along the equatorial axis (Fig. 3B) was almost insignificant, we observed a considerable difference when comparing the upstream and downstream concentration profiles (Fig. 3, D and E). For any Reynolds numbers in the Re range examined here, both the O2 concentrations and their gradients downstream were smaller than those upstream. The effect of increasing flow on the O2 distribution along the z axis was found to be similar to that along the x axis, that is, O2 concentration increased with an increase of Reynolds number.

An important phenomenon, previously poorly known, appeared (Fig. 4, A and B) when plotting the O2 flux and O2 concentration at the aggregate’s surface as a function of angle Θ with respect to the x axis (see Fig. 1A for the definition of the angle Θ; we note that the angle Θ does not follow the standard definition in the spherical coordinate system). Here, the term flux refers to J = uCDC. However, since we consider impermeable aggregates in this study, we imposed the no-slip boundary condition (u = 0) at the aggregate’s surface and reduced the interfacial flux to the diffusive flux, −DC. The flux was calculated using the standard lattice Boltzmann (LB) procedure [see (15) for the further details]. As demonstrated here, one can consider two regions: (i) the lower half of the aggregate (Θ < 0°) where both the flux and the concentration remain almost unchanged and (ii) the upper half of the aggregate where O2 concentration profiles decrease near the equator and reach a minimum at Θ = +90°. The simultaneous variation of concentration and its gradient (flux) with Θ demonstrates that the biological activity within a porous aggregate can be modeled adequately neither by assuming constant concentration (Dirichlet) condition nor by assuming constant flux (Neumann) condition at the aggregate’s surface. This was also recently found in an experimental study measuring O2 flux in the downstream and equatorial region of yeast-agar spheres used as model aggregates (14). Here, we would like to highlight an observed nonintuitive trend: While the O2 concentration increased at the downstream pole (Θ = +90°) with increasing Re, the corresponding flux decreased (Fig. 4, A and B). This phenomenon was associated with the fact that, at the downstream region, the flow direction (and thus advective flux) opposed the direction of diffusive flux. In Fig. 4C, we have plotted the effective diffusive boundary layer thickness, δ, as a function of Θ (see Fig. 6A and Materials and Methods for the definition of δ). As shown, δ decreased with increasing Reynolds numbers in the entire range of Θ. However, this effect was mitigated in the downstream region, particularly at Θ = +90°. This explains why the diffusive flux downstream is lower than the equatorial and upstream flux.

Fig. 4 Flux, concentration and diffusive boundary layer thickness at the aggregate’s surface for various Re.

O2 flux (A), O2 concentration (B), and effective diffusive boundary layer thickness (C) at the aggregate’s interface as a function of Θ for different Reynolds numbers. The color code holds for all panels.

At the steady state, the actual total O2 uptake by all heterotrophic organisms within an aggregate, Qact, is supplied by the incoming O2 flux across the entire surface of the aggregate. This implies that the volume integral of the reaction term is equivalent to the surface integral of flux into the aggregate, yieldingEmbedded Image(1)

Since the angle between the vector J and the surface element ds is not known a priori, it is more convenient to calculate the second integral for computation of Qact. Accordingly, one may define a normalized average flux Javr = (Qact/S) where S is the surface area of the aggregate. We computed the total O2 uptake and the average O2 flux Javr for the entire range of Re ≤ 10 (Fig. 5A). A close inspection of the computed values reveals that the enhancement of the total O2 uptake as a function of Reynolds number was limited to less than 1% in our study, suggesting a reaction-limited process. Next, we quantified the relative errors of equatorial flux, |JxJavr |/Javr, and downstream flux, |JzJavr |/Javr, with respect to the average O2 flux, as shown in Fig. 5B. The error associated with the total O2 uptake amounts to ca. 7% when the total surface area is multiplied by Jx. Note that the error percentages reported here depend on C and on the associated Km and Vm (10, 14). Although the equatorial flux, Jx can be considered as a relatively good approximation of Javr, measuring the equatorial O2 concentration profiles is rather a highly challenging task. Downstream measurements of O2 concentration profiles are easier to perform but lead to underestimation of Javr by up to 35% (Fig. 5B). For example, for an aggregate sinking with a typical Reynolds number of Re = 2, the calculation of the total O2 uptake based on downstream measurement of O2 concentration profiles results in almost 32% underestimation of the actual value. Except in a recent experiment (14) on perfect spheres with yeast as model aggregates, this fact was unknown.

Fig. 5 Effect of sinking velocity on the total uptake, downstream and equatorial fluxes and comparison of simulation results with theoretical prediction for the uptake enhancement.

(A) Total O2 uptake (red color) and the average flux Javr (blue color) versus Reynolds number. (B) The relative error of Jx (red color) and Jz (blue color) with respect to Javr. (C) The ratio of the flux measured at downstream to the average flux (Ω, red color) and the fitted model (Eq. 2, blue color) as a function of Re. (D) Uptake enhancement due to the advection. The red curve is Eq. 4, and the blue line is the calculated total uptake using the present model.

It is also interesting to look at the ratio of downstream to the average flux, Ω = Jz/Javr, as a function of Reynolds number shown in Fig. 5C. In the absence of advection, Ω is equal to unity and runs into an asymptotic behavior for large Re. The corresponding data are best represented by a function of the formEmbedded Image(2)where A = 0.6878 and B = −0.05715, as plotted in the Fig. 5C. When considering the range of 1 ≤ Re ≤ 10, Ω varied slightly from 0.71 to 0.65, which is close to the value of A. On the basis of this finding, one can obtain the actual (correct) O2 uptake asEmbedded Image(3)when the downstream flux is used.

The uptake enhancement due to advection can be presented by Sherwood number Sh as the ratio of total mass transport into/out of a moving aggregate (advection + diffusion + reaction) to that of a nonmoving aggregate (diffusion + reaction). For a solid aggregate model with the constant concentration condition imposed at the aggregate’s surface, the O2 uptake enhancement follows the relationEmbedded Image(4)where Sc = ν/D is the Schmidt number (11). The comparison of the uptake enhancement for a solid aggregate model versus the porous aggregate of the present study revealed that imposing constant concentration at the aggregate’s surface overestimated the total O2 uptake by 1.8- to 15-fold within the Reynolds number range of 0.01 ≤ Re ≤ 10 (Fig. 5D).

The application range of the model presented here is not limited to the investigation of the O2 consumption of sinking aggregates only. This model can address a variety of microprocesses in the microenvironment of sinking porous aggregates and phytoplankton colonies, which are discussed below.


Small-scale fluxes and the microenvironment of nutrients and gases (for example, NH4+, CO2, and O2) in phytoplankton (16, 17), phytoplankton colonies (18, 19), and sinking particles and aggregates (10, 11, 20) have been in focus for decades. Pure theoretical studies of small-scale fluxes to phytoplankton, phytoplankton colonies, and small particles have often assumed Dirichlet boundary conditions (1620). In contrast, empirical studies of small-scale O2 fluxes at the aggregate-water interface measured by O2 microsensors have often been analyzed and compared with mass transfer theory assuming Neumann boundary conditions (8, 10, 11, 14, 21, 22). The present study shows that neither of these boundary conditions accurately describes small-scale fluxes to and from aggregates at Re > 0.1. Furthermore, our mechanistic understanding of small-scale fluxes and pelagic microenvironments has been constrained by the lack of a model that accurately describes mass transfer and reaction within sinking porous aggregates incorporating Michaelis-Menten kinetics of the associated biota as variables. The model presented here includes these variables and accurately reproduced measured O2 fluxes and O2 distributions associated to sinking porous aggregates (10, 14). Thus, this model represents a novel tool by which we can explore physical, chemical, and biological constraints of pelagic microenvironments. Below, we present a few examples of emergent research questions in plankton ecology, microbial ecology, and biogeochemistry to which the model has a substantial potential to contribute with significant mechanistic understanding in near future studies.

Phytoplankton colonies are pelagic microenvironments with steep gradients of O2, pH, and nutrients due to their CO2 fixation and cellular growth. Consequently, the growth conditions for associated biota are significantly different from those in the ambient water. Trichodesmium is a significant N2-fixing organism, which often occurs as colonies in the (sub-)tropical ocean where it is a key organism in N2 fixation and CO2 sequestration from the atmoshpere (23). An early study of the O2 microenvironment in Trichodesmium colonies suggested that these were generally anoxic in their center (18). Later studies, however, have reported that O2 concentrations are, in contrast, supersaturated within colonies during photosynthesis and that the O2 concentration remains higher than 50% air saturation during darkness (19, 24). The nitrogenase enzyme is highly sensitive to O2 (25), and the ability of Trichodesmium to do N2 fixation despite of high O2 concentrations within colonies is consequently challenging (19). The recent introduction of nanoscale secondary ion mass spectrometry (nanoSIMS) has enabled detailed studies of cellular growth at a single-cell level. These nanoSIMS studies of Trichodesmium colonies have revealed a high diversity of cell-specific N2 fixation and C fixation at a single-cell level (19, 23), which may be explained by gradients of O2 and pH within colonies, as demonstrated with empirical measurements (19). Thus, microsensor analysis, in combination with nanoSIMS, and modeling represent a promising new approach to revealing the constraints that determine actual cellular growth rates and microscale growth conditions inside N2-fixing cyanobacterial colonies (19, 26). N2-fixing cyanobacterial colonies have been reported to release newly fixed N2 fixation as NH4+ and dissolved organic nitrogen (26, 27)—a loss that may appear contradictory to the fact that N2 fixation is a highly energy-demanding process. It has been suggested that coexisting heterotrophic bacteria gaining NH4+ from the cyanobacteria within colonies may modulate the chemical microenvironment of both O2 and pH through their respiration to improve growth conditions of the cyanobacteria including N2 fixation and trace metal uptake (28). Unfortunately, ammonium microsensors are not available for measurements in marine environments because of interference with salt (29). Instead, NH4+ fluxes from colonies and aggregates have been measured in the ambient water, and the concentration distributions have been modeled on the basis of the measured fluxes and aggregate size and on assumptions of diffusivity (14, 26). This model demonstrated that NH4+ concentration may be <60-fold higher inside these colonies, as compared to concentrations in the ambient water (26). The NH4+ microenvironment in N2-fixing cyanobacterial colonies supports not only biota living within the colony interior but also other N-limited phytoplankton living in the ambient water (30). Using stable isotope tracers, it has been demonstrated that N2 fixation and the small-scale NH4+ fluxes from cyanobacterial colonies to the ambient water support both the microbial and the classical food web in the sea (30). The newly fixed N2 even supports the growth of chain-forming diatoms (for example, Chaetoceros), which form fast-sinking and transport carbon and nutrients to depth in the ocean (30, 31). Obviously, the chemical microenvironment of Trichodesmium is of great importance for both small-scale N2 fixation and large-scale carbon sequestration in the ocean. Nevertheless, the growth conditions and their constraints within these colonies remain enigmatic (19) but may be understood through modeling followed by design of experiments in the future.

Recently, the microbial communities and potential denitrification activities associated with Trichodesmium have attracted much attention in the research community (32, 33). The presence of potentially denitrifying microorganisms associated with Trichodesmium has raised the hypothesis that Trichodesmium may even support denitrification in its O2 microenvironment (34), as recently demonstrated experimentally in N2-fixing Nodularia colonies (21). In the latter study, it was shown that a complete aerobic and anaerobic N cycle can be present in the microenvironment of Nodularia colonies. Denitrification and dissimilatory nitrate reduction to ammonium within colonies were limited by nitrification within the colony or small-scale nitrate fluxes from the ambient water to the anoxic colony center (21). Low nitrification activities have also been observed in diatom aggregates (14, 22), but its reasons remain poorly understood. Export production, including sinking aggregates, connects aerobic processes in the upper ocean with anaerobic processes in the hypoxic oxygen minimum zones (35) with ambient O2 concentrations <20 μM where most aggregates and fecal pellets are predicted to be anoxic (10, 14, 21, 22). Recently, it was demonstrated that sinking (diatom) aggregates and fecal pellets can be important sites connecting aerobic N-transformation processes in the euphotic zone with the anaerobic N cycle in hypoxic oxygen minimum zones where ambient nitrate concentrations are high and N losses due to N2 production through denitrification and anaerobic ammonium oxidation (anammox) are high (10, 14, 21, 22). Using the model of the present study, we can now examine physical, chemical, and biological constraints of microbial coexistence and the coupling of anoxic microenvironments and the aerobic/anaerobic N cycle in the ocean.

Sinking aggregates link small-scale processes with large-scale biogeochemical fluxes in the ocean. Their role in the biological carbon pump is largely governed by physical aggregation and disaggregation processes, their microbial transformation, and their role as food source for zooplankton during descent (36). Intriguingly, because of the hydrodynamics and small-scale fluxes, as described by the present model, the microbial hydrolysis and respiration processes leave a plume of DOM and nutrients behind sinking aggregates, which is detected by motile, colonizing bacteria and by zooplankton (for example, copepods) (11, 37, 38). Thus, microbial colonization and transformation and zooplankton grazing interact and can be tightly coupled on sinking marine aggregates. The constraints of these co-occurring processes include the aggregate’s chemical composition as nutritional value, its porosity and permeability, sinking velocity, and their interactions and biological activities (for example, hydrolysis and respiration) (7), while the concentrations and distribution of DOM plume must be largely governed by small-scale diffusion and advection at the rear of the aggregate-water interface. These parameters can all be included in the present model. The combination of modeling with empirical measurements of activities using stable and radioactive isotope tracers, microsensors in combination with nanoSIMS, and CardFISH (Catalyzed reporter deposition Fluorescence In Situ Hybridization) represents an exciting and promising novel approach for quantitative and mechanistic understanding of these multiple parameters and their complexity coupling small-scale and large scale fluxes in the ocean.


Flow field

To model the flow field around a sinking porous aggregate, we considered an impermeable sphere of radius R sinking with a constant velocity U0. The steady downward motion of an aggregate in the stagnant ambient water is equivalent to the flow past a stagnant aggregate subjected to a constant upward velocity U in the far field (UU0). In this study, we considered the latter case (Fig. 1A). Because of the axial symmetry around the z axis, one may consider the flow and concentration fields as a two-dimensional problem where the coordinate system has its origin at the aggregate’s center. The governing equation for the flow field around the aggregate is given by the incompressible Navier-Stokes equationEmbedded Image(5)where u, t, ρ, ν, and p are fluid velocity vector, time, fluid density, kinematic viscosity, and fluid pressure, respectively. In this study R, ν, and ρ were fixed to 1.75 mm, 1.16 × 10−2 cm2/s, and 1.025 g/cm3, respectively.

Concentration field

The governing advection-diffusion reaction equation for solute concentration within the considered domain reads asEmbedded Image(6)

In the above equation, ε, C, and D represent aggregate porosity, concentration, and diffusion coefficient of the solute (or nutrient) in the ambient water, respectively. The reaction term Rr is associated with consumption (or production) of solutes via microbial activities and is modeled via Michaelis-Menten kinetics. In the reaction term, Vm denotes the maximum reaction rate achieved at saturating solute concentration, and Km is the half-saturation constant. Inclusion of the reaction term allows us to address reaction-limited versus diffusion-limited processes. Solute concentration sufficiently far from the aggregate is set to the ambient concentration C. The maximum reaction rate Vm can be more conveniently represented using the Damköhler number. Since we focus on impermeable aggregates, the appropriate Damköhler number Da is defined as the ratio of biological reaction rate to the diffusive mass transfer rate. This implies that Da increases either because of a decrease in the diffusive mass transfer rate or because of an increase in the biological reaction rate and can be equivalently written as Da = τDR where τD and τR denote diffusion and reaction time scales, respectively. The diffusion time scale is given by R2D−1; thus, the maximum reaction rate can be estimated as VmDaCDR−2. Hereby, the biological activities inside an aggregate can be characterized by Da and Km.

Numerical method

To solve the above equations, the LB method (15, 39) was used, and extensive computer simulations were performed. The core of the current LB method was explained and verified previously (13). The current model is an extension of the previous model (13) by inclusion of the reaction term in the concentration equation (Eq. 6).

Model validation

Results of the LB model were previously verified (13). To verify the additional reaction term, the effective diffusive boundary layer thickness δ is shown in Fig. 6A for the case of Re = 0. The effective diffusive boundary layer is an effective thickness, defined as the distance from the aggregate’s surface to the point where the tangent line to the concentration profile at the aggregate’s surface crosses the ambient concentration C. This way, the diffusive flux is given by D(CC0)/δ, with C0 as the concentration at the aggregate’s surface. As demonstrated, for the case of Re = 0, the δ for the O2 concentration is practically equal to the radius of the aggregate. One can theoretically obtain the same result (δ/R = 1) from Crank’s equation for solute concentration around a sphere (8, 40). A further verification is provided by comparing our results with experimental measurements of concentration field within and around a sinking aggregate (10), as shown in Fig. 6 (B and C, respectively). Although the shape of the marine aggregates investigated by Ploug (10) deviated from perfect spheres (Fig. 6B), the present model was able to reproduce similar pattern and quantity of the O2 concentration isolines in the central xz plane (Fig. 6C). It is noted that, in contrast to the modeled aggregate, the isolines in the experiment are a bit shifted to the rear of the aggregate. This difference might be due to the differences between the shape of the modeled and the real aggregate.

Fig. 6 Calculation of the diffusive boundary layer thickness for the case of Re = 0 and comparison of measured and simulated O2 concentration fields.

(A) The tangent line to the concentration profile at the interface crosses the ambient concentration C at a distance of R from the interface, that is, δ/R = 1, for the case of Re = 0. (B) The concentration isolines of a sinking porous aggregate (with a sinking velocity of 64 m/day) measured in the experiment [the figure was taken from a previous study of Ploug (10)]. (C) The concentration isolines of a sinking porous aggregate, simulated using the current model with Km = 1 and Da = 1.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: A.K. received a partial funding from Max Planck Society (MPG) and a grant from Deutsche Deutsche Forschungsgemeinschaft (DFG) (KH 31/8-1). H.P. received funding from University of Gothenburg and grants from the Swedish Research Council (VR, Dnr: 621-2011-4406 and Dnr: 2015-05322). M.I. received funding from the Helmholtz Young Investigator Group SEAPUMP “Seasonal and regional food web interactions with the biological pump” (VH-NG-1000), the Alfred Wegener Institute for Polar and Marine Research, and the DFG-Research Center/Cluster of Excellence MARUM “The Ocean in the Earth System.” We thank both anonymous referees for valuable comments. Author contributions: N.M. extended, parameterized, and optimized the computer program, performed the simulations, provided and analyzed the results, and wrote the manuscript. B.L. provided the first version of the computer program. M.I. contributed to the link to the biology and wrote the manuscript. M.M.K. reviewed the manuscript and made valuable comments. H.P. provided the link between the numerical results and biology and wrote the manuscript. A.K. conceived the original idea of this study, led the entire study, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article