## Abstract

N_{2}-fixing colonies of cyanobacteria and aggregates of phytoplankton and detritus sinking hundreds of meters per day are instrumental for the ocean’s sequestration of CO_{2} 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 O_{2} 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 O_{2} flux in these aggregates by as high as 15-fold.

## INTRODUCTION

The rising CO_{2} concentration in the Earth’s atmosphere is partly counterbalanced by CO_{2} sequestration in phytoplankton and N_{2}-fixing cyanobacteria. One of the major topics in aquatic microbial ecology deals with the fate of fixed CO_{2} 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 (*4*–*6*), 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 (*7*–*9*). 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 O_{2} uptake by individual aggregates are usually made either by laboratory experiments or by mathematical models based on mass transfer theories. Experimental determination of O_{2} uptake is often based on calculation of O_{2} flux using O_{2} 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 the aggregate’s surface area to obtain the total uptake (*2*, *8*, *10*). In an attempt to quantify the O_{2} 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 O_{2} 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 O_{2} field within and around sinking marine aggregates and accurately quantify total O_{2} uptake by microbes associated with the aggregates.

## RESULTS

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 *U*_{0}, the corresponding Reynolds number reads *Re* = *U*_{0}*R*/υ, 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 (*U*_{∞} ≡ *U*_{0}) 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 *Z*_{U10}, as a function of *Re* (Fig. 1D).

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 O_{2} diffusion coefficient, ambient concentration, porosity, half-saturation constant, and Damköhler number to *D* = 1.70 × 10^{−5} cm^{2}/s, *C*_{∞} = 300 μM, ε = 0.95, *K*_{m} = 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 *K*_{m} = 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 O_{2} 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 O_{2} 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 *V*_{m} is proportional to *Da*. For the case of *Da* = 10, the microbial O_{2}, respiration rate is much larger than the diffusive rate, and thus, the O_{2} uptake by the aggregate is diffusion limited (Fig.1E). Furthermore, it is interesting to examine the effect of *K*_{m} and ε variation on the O_{2} concentration profile when other parameters remain unchanged. Although the O_{2} concentration increases with the increase of *K*_{m}, the concentration variation is not significant as long as *K*_{m} is sufficiently lower than *C*_{∞} (Fig. 1F). Since the interstitial water inside the aggregate is reduced with the decrease of porosity, the O_{2} concentration inside the aggregate decreases too (Fig. 1G). Figure 2 illustrates the normalized O_{2} 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 O_{2} concentration of *C*/*C*_{∞} = 1.0 falls within a distance of approximately 3*R* 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 O_{2} concentration profiles along the *x* (horizontal) and *z* (vertical) axes is shown in Fig. 3. Those computations revealed that the O_{2} 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, O_{2} 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 O_{2} concentrations closely approached the ambient concentration decreased with increasing flow. The O_{2} concentration profiles, however, become more compressed with increasing *Re* (for example, compare *Re* = 5 and *Re* = 10).

Unlike O_{2} concentration along the *x* axis, concentration along the *z* axis was asymmetric (Fig. 3C). While the difference between the upstream O_{2} 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 O_{2} concentrations and their gradients downstream were smaller than those upstream. The effect of increasing flow on the O_{2} distribution along the *z* axis was found to be similar to that along the *x* axis, that is, O_{2} concentration increased with an increase of Reynolds number.

An important phenomenon, previously poorly known, appeared (Fig. 4, A and B) when plotting the O_{2} flux and O_{2} 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** =

**u**

*C*+

*D*∇

*C*. 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,

*D*∇

*C*. 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 O

_{2}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 O

_{2}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 O

_{2}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.

At the steady state, the actual total O_{2} uptake by all heterotrophic organisms within an aggregate, *Q*_{act}, is supplied by the incoming O_{2} 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, yielding(1)

Since the angle between the vector ** J** and the surface element

*d*

**is not known a priori, it is more convenient to calculate the second integral for computation of**

*s**Q*

_{act}. Accordingly, one may define a normalized average flux

*J*

_{avr}= (

*Q*

_{act}/

*S*) where

*S*is the surface area of the aggregate. We computed the total O

_{2}uptake and the average O

_{2}flux

*J*

_{avr}for the entire range of

*Re*≤ 10 (Fig. 5A). A close inspection of the computed values reveals that the enhancement of the total O

_{2}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, |

*J*

_{x}−

*J*

_{avr}|/

*J*

_{avr}, and downstream flux, |

*J*

_{z}−

*J*

_{avr}|/

*J*

_{avr}, with respect to the average O

_{2}flux, as shown in Fig. 5B. The error associated with the total O

_{2}uptake amounts to ca. 7% when the total surface area is multiplied by

*J*

_{x}. Note that the error percentages reported here depend on

*C*

_{∞}and on the associated

*K*

_{m}and

*V*

_{m}(

*10*,

*14*). Although the equatorial flux,

*J*

_{x}can be considered as a relatively good approximation of

*J*

_{avr}, measuring the equatorial O

_{2}concentration profiles is rather a highly challenging task. Downstream measurements of O

_{2}concentration profiles are easier to perform but lead to underestimation of

*J*

_{avr}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 O

_{2}uptake based on downstream measurement of O

_{2}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.

It is also interesting to look at the ratio of downstream to the average flux, Ω = *J*_{z}/*J*_{avr}, 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 form(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) O_{2} uptake as(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 O_{2} uptake enhancement follows the relation(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 O_{2} 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 O_{2} 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.

## DISCUSSION

Small-scale fluxes and the microenvironment of nutrients and gases (for example, NH_{4}^{+}, CO_{2}, and O_{2}) 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 (*16*–*20*). In contrast, empirical studies of small-scale O_{2} fluxes at the aggregate-water interface measured by O_{2} 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 O_{2} fluxes and O_{2} 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 O_{2}, pH, and nutrients due to their CO_{2} fixation and cellular growth. Consequently, the growth conditions for associated biota are significantly different from those in the ambient water. *Trichodesmium* is a significant N_{2}-fixing organism, which often occurs as colonies in the (sub-)tropical ocean where it is a key organism in N_{2} fixation and CO_{2} sequestration from the atmoshpere (*23*). An early study of the O_{2} microenvironment in *Trichodesmium* colonies suggested that these were generally anoxic in their center (*18*). Later studies, however, have reported that O_{2} concentrations are, in contrast, supersaturated within colonies during photosynthesis and that the O_{2} concentration remains higher than 50% air saturation during darkness (*19*, *24*). The nitrogenase enzyme is highly sensitive to O_{2} (*25*), and the ability of *Trichodesmium* to do N_{2} fixation despite of high O_{2} 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 N_{2} fixation and C fixation at a single-cell level (*19*, *23*), which may be explained by gradients of O_{2} 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 N_{2}-fixing cyanobacterial colonies (*19*, *26*). N_{2}-fixing cyanobacterial colonies have been reported to release newly fixed N_{2} fixation as NH_{4}^{+} and dissolved organic nitrogen (*26*, *27*)—a loss that may appear contradictory to the fact that N_{2} fixation is a highly energy-demanding process. It has been suggested that coexisting heterotrophic bacteria gaining NH_{4}^{+} from the cyanobacteria within colonies may modulate the chemical microenvironment of both O_{2} and pH through their respiration to improve growth conditions of the cyanobacteria including N_{2} fixation and trace metal uptake (*28*). Unfortunately, ammonium microsensors are not available for measurements in marine environments because of interference with salt (*29*). Instead, NH_{4}^{+} 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 NH_{4}^{+} concentration may be <60-fold higher inside these colonies, as compared to concentrations in the ambient water (*26*). The NH_{4}^{+} microenvironment in N_{2}-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 N_{2} fixation and the small-scale NH_{4}^{+} fluxes from cyanobacterial colonies to the ambient water support both the microbial and the classical food web in the sea (*30*). The newly fixed N_{2} 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 N_{2} 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 O_{2} microenvironment (*34*), as recently demonstrated experimentally in N_{2}-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 O_{2} 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 N_{2} 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.

## MATERIALS AND METHODS

### 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 *U*_{0}. 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 (*U*_{∞} ≡ *U*_{0}). 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 equation(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} cm^{2}/s, and 1.025 g/cm^{3}, respectively.

### Concentration field

The governing advection-diffusion reaction equation for solute concentration within the considered domain reads as(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 *R*_{r} is associated with consumption (or production) of solutes via microbial activities and is modeled via Michaelis-Menten kinetics. In the reaction term, *V*_{m} denotes the maximum reaction rate achieved at saturating solute concentration, and *K*_{m} 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 *V*_{m} 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* = τ_{D}/τ_{R} where τ_{D} and τ_{R} denote diffusion and reaction time scales, respectively. The diffusion time scale is given by *R*^{2}*D*^{−1}; thus, the maximum reaction rate can be estimated as *V*_{m} ≈ *DaC*_{∞}*DR*^{−2}. Hereby, the biological activities inside an aggregate can be characterized by *Da* and *K*_{m}.

### 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*(*C*_{∞} − *C*_{0})/δ, with *C*_{0} as the concentration at the aggregate’s surface. As demonstrated, for the case of *Re* = 0, the δ for the O_{2} 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 O_{2} 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.

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.

## REFERENCES AND NOTES

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

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).