Theoretical Insight into the Biodegradation of Solitary Oil Microdroplets Moving through a Water Column

In the aftermath of oil spills in the sea, clouds of droplets drift into the seawater column and are carried away by sea currents. The fate of the drifting droplets is determined by natural attenuation processes, mainly dissolution into the seawater and biodegradation by oil-degrading microbial communities. Specifically, microbes have developed three fundamental strategies for accessing and assimilating oily substrates. Depending on their affinity for the oily phase and ability to proliferate in multicellular structures, microbes might either attach to the oil surface and directly uptake compounds from the oily phase, or grow suspended in the aqueous phase consuming solubilized oil, or form three-dimensional biofilms over the oil–water interface. In this work, a compound particle model that accounts for all three microbial strategies is developed for the biodegradation of solitary oil microdroplets moving through a water column. Under a set of educated hypotheses, the hydrodynamics and solute transport problems are amenable to analytical solutions and a closed-form correlation is established for the overall dissolution rate as a function of the Thiele modulus, the Biot number and other key parameters. Moreover, two coupled ordinary differential equations are formulated for the evolution of the particle size and used to investigate the impact of the dissolution and biodegradation processes on the droplet shrinking rate.


Introduction
After a natural or accidental release of crude oil in the sea, part of the oil ends up in the form of droplets moving through the seawater column. The droplets may be created either at the sea surface during the breakup of an oil slick (i.e., floating oil layer) by sea waves [1,2], or at the seafloor during the atomization of live crude oil (i.e., gas/oil mixture) extruding at sufficiently high speed from a natural crack or a broken wellhead [3][4][5]. The latter case occurred, for example, after the blowout of the Deepwater Horizon rig in the Gulf of Mexico where the addition of the chemical dispersant Corexit in the leaking crude oil resulted in clouds of droplets travelling underwater along with sea currents [6,7]. At present, there are no practical means for the collection or in situ treatment of oil droplets in vast bodies of marine waters and, inevitably, their removal relies solely on natural attenuation processes, notably on dissolution and biodegradation. Specifically, it is anticipated that in the long run, most of the released oil in the sea is consumed by autochthonous oil-degrading microorganisms (bacteria, fungi, yeasts) that have developed appropriate machinery for accessing and assimilating oily substrates [8][9][10]. In this way, crude oil enters as a nutrient into the marine food for the fate of oil droplets in marine waters account only for the direct uptake strategy [32][33][34], neglecting any effects of biodegradation in the bulk aqueous phase or the formation of biofilm over the droplet and bioreaction therein.
In this work, a compound particle model (CPM) is developed for the biodegradation of solitary oil microdroplets moving through a water column. The compound particle is of the core-shell type and consists of an oily core that is successively surrounded by a bioreactive skin of negligible thickness and another bioreactive shell of finite thickness (Figure 1). The bioreactive skin represents a thin layer of microbes that uptake oil directly from the oily phase, whereas the bioreactive shell represents a distinct biofilm phase. In line with the abovementioned microbial strategies of oil uptake, the model accounts for all three modes of biodegradation: direct interfacial uptake, bioreaction in the bulk aqueous phase, and bioreaction in a biofilm formed around the droplet. A set of simplifying hypotheses is introduced so as to make the mathematical analysis tractable, and the governing equations solvable by analytical methods. The most important hypotheses are that the compound particle is considered to move as a non-deforming rigid sphere, the flow of the aqueous phase is dominated by viscous stresses, and the transport of dissolved oil in the biofilm phase is dominated by diffusion, whereas in the bulk aqueous phase it is dominated by advection. The analysis of the local mass balances results in a closed-form expression for the overall dissolution rate as a function of the Biot number, the Thiele modulus, the thickness of the biofilm, and the diffusivity and solubility ratios. Furthermore, from the overall mass balances, two coupled ordinary differential equations are established for the evolution of the particle size.

Model Formulation
With reference to Figure 1, the process under consideration is the transport and reaction of dissolved oil, denoted as the A solute, from the oil droplet (Ωλ) to the surrounding biofilm (Ωβ) and aqueous (Ωυ) phases. The thick line at the oil-biofilm interface (Sλβ) represents a thin layer of microbes that uptake oil compounds directly from the oily phase. The first step in the theoretical analysis is to determine the oil dissolution rate at the droplet surface, based on an appropriate formulation of the local mass balances (Section 2.3). The second step is to determine the droplet shrinking rate using the overall mass balances (Section 2.4). Before proceeding with the mathematical

Model Formulation
With reference to Figure 1, the process under consideration is the transport and reaction of dissolved oil, denoted as the A solute, from the oil droplet (Ω λ ) to the surrounding biofilm (Ω β ) and aqueous (Ω υ ) phases. The thick line at the oil-biofilm interface (S λβ ) represents a thin layer of microbes that uptake oil compounds directly from the oily phase. The first step in the theoretical analysis is to determine the oil dissolution rate at the droplet surface, based on an appropriate formulation of the local mass balances (Section 2.3). The second step is to determine the droplet shrinking rate using the overall mass balances (Section 2.4). Before proceeding with the mathematical analysis, certain key considerations on modeling the different biodegradation modes (Section 2.1) and a set of basic hypotheses (Section 2.2) are set forth.

Considerations on Modeling the Three Major Biodegradation Modes
A few remarks are in order with regard to the theoretical modeling of each one of the three basic modes of biodegradation; that is, direct interfacial uptake, bioreaction in the bulk aqueous phase, and bioreaction in a biofilm formed around the droplet.
The first type of oil-degrading microbes is the flatlanders; that is superhydrophobic microbes able to firmly adhere to the oil surface and directly uptake organic compounds from the oily phase. Here, it is considered that the oil surface (S λβ ) is fully and uniformly covered by flatlanders. Partial coverage is expected to lead to more complex phenomena of fluid dynamics and solute transport and, thus, deserves to be investigated separately. The layer of flatlanders is usually found embedded in the oil side of the oil-water interface [35] and can be viewed as a bioreactive skin (interphase) of negligible thickness (~a few µm) on the droplet scale of observation (~100 µm). As the microbes have direct access to the oily substrate, the oil consumption rate is considered to be limited only by the intrinsic microbial kinetics. Under these conditions, this mode of biodegradation is essentially decoupled from the dissolution of oil to the surrounding phases. The physical presence of microbes on the droplet surface and the process of interfacial reaction are assumed to affect only implicitly the dissolution of oil; that is, by (possibly) changing the value of oil solubility.
The second type of oil-degrading microbes is the drifters; that is, hydrophilic microbes that remain suspended in the bulk aqueous phase (Ω υ ) and consume solubilized (molecular or micellar) oil. For this biodegradation mode, it is considered that the concentration of microbes is constant throughout the aqueous phase and the oil consumption rate follows first-order kinetics. In addition, solute A represents both molecular and micellar oil and, thus, the action of surfactants is taken into account only implicitly by modifying the apparent solubility of oil in the aqueous phase.
The third type of oil-degrading microbes is the biofilm formers; that is, microbes able to actively construct three-dimensional biofilm communities over the oil surface. The thickness of the biofilm might be appreciable and, thus, the biofilm is viewed as a distinct phase on the droplet scale of observation. Supplementary hypotheses for this mode include a uniform biofilm thickness, constant concentration of active microbes within the biofilm, and first order kinetics for the oil consumption rate. Interstitial flow is neglected and solute transport within the biofilm is dominated by diffusion.
With regard to the microbial proliferation rate, it is customary to assume a linear dependence on the concentration of active cells, that is where α denotes the physical domain in which the microbes grow and takes the values α = β, υ, λβ. All of the primary symbols are defined in the nomenclature. The tilde (~) over a variable or parameter denotes a dimensional quantity, whereas the lack of it denotes a dimless quantity. The term dimensionless is abbreviated to dimless throughout the paper. The specific growth rate µ C,α is usually considered to follow Monod kinetics.

of 29
Any possible effects of lag phase, cell maintenance, limitation by electron acceptors, nitrogen and phosphorous sources, substrate cometabolism or inhibition are neglected. Of particular interest are the limiting forms for the specific growth rate under sufficiently high or low concentration.
The zeroth-order kinetics is expected to be applicable in the interfacial uptake mode because the microbial cells have access to the pure oily phase. On the other hand, the first-order kinetics is expected to be applicable in the suspended and biofilm growth modes because of the low solubility of oil compounds in aqueous phases. In all cases, the oil consumption rate is considered to be proportional to the cell proliferation rate. Therefore, the volumetric consumption rate of dissolved oil in a bulk phase, is given by for α = β, υ; and the surficial consumption rate on the droplet surface, is given by Here, B α and B λβ is the volumetric and surface concentration of cells, respectively. The minus sign denotes consumption rates. In this work, the cell concentration is treated as a constant for all three biodegradation modes.

Basic Hypotheses for the Hydrodynamics and Mass Transport
In addition to the previous considerations for the biodegradation process, a set of hypotheses is introduced for the flow and mass transport processes so as to simplify the mathematical description as much as possible while retaining the most important mechanisms. First, the external flow in the unbounded aqueous phase is dominated by viscous stresses and, thus, characterized by a low Reynolds number (Re υ = R P U ρ υ / µ υ 1). Second, the internal recirculating flow and the deformation of the particle are considered to be negligible. In all cases, the adsorption of biopolymers and microbial cells to the oil-water interface is expected to hinder the interfacial mobility and, consequently, diminish the internal flow in the oily phase. On the basis of a combination of small particle size, slow velocity and rigid-like interface, it is expected that interfacial tension dominates over viscous and gravitational forces that tend to deform the particle and the system is, thus, characterized by low capillary and Bond numbers (Ca = µ υ U/ γ βυ 1, Bo = R 2 P ∆ ρ g/ γ βυ 1; where ∆ ρ = ρ υ − ρ p is the excess density and γ βυ is the interfacial tension at the particle surface) [36]. Therefore, the particle, either simple or compound, is considered to move as a rigid sphere. Third, the transport of dissolved oil in the biofilm phase is dominated by diffusion and, thus, characterized by a low Péclet number (Pe β = R P U β / D Aβ 1). On the other hand, mass transport in the bulk aqueous phase is considered to be dominated by advection and characterized by a high Péclet number (Pe υ = R P U/ D Aυ 1). In both phases, solute diffusion is considered to obey Fick's constitutive law. Fourth, the oily phase is treated as a single compound and mass transfer therein is not taken into account (e.g., the solute A represents the total petroleum hydrocarbon in the case of crude oil). Finally, the quasi-steady state hypothesis is adopted for the analysis of the flow and mass transport problems at the local level (Section 2.3). Besides a high Péclet and a low Reynolds number, this assumption also requires a low droplet shrinking rate as compared to the characteristic velocity of the external flow. Thereafter, the evolution of the particle size is treated as a sequence of steady states in Section 2.4.

Overall Dissolution Rate: Analysis of the Local Mass Balances
Under the detailed set of considerations and hypotheses given in the previous subsections, mass transport is described in the context of the CPM by the following equations It is possible to further reduce the complexity of the above set of governing equations by introducing two educated hypotheses. First, the tangential diffusion in the spherical shell is neglected. Strictly, this hypothesis holds for a thin shell ( δ β R P ) or fast reaction (Da β 1). Thus, Equation (6b) becomes Second, the continuity of the mass flux at the υβ-interface is imposed in an average sense by demanding equality of the surface-averaged fluxes, instead of equality of the local fluxes. Therefore, Equation (6c) is expressed as follows The above equation can be tidied up by considering that S βυ is a spherical surface at r = R P with d S = r 2 sin θdθdϕ and n βυ = e r . Also, the radial mass flux in the β-phase is independent of the polar and azimuthal angles, while the surface averaged flux in the right hand side of Equation (8a) defines the dissolution rate from the particle surface to the υ-phase. Thus, Equation (8a) becomes The partition coefficient of oil at the υβ-interface, H A,υ/β , is approximated as the solubility ratio in the corresponding phases. By replacing Equations (6b) and (6c) with Equations (7) and (8c), respectively, and also by introducing dimless quantities, the mass transport problem defined in Equations (6a)-(6f) obtains the form c Aβ (R C ) = 1, at r = R C , (9a) For the non-dimensionalization, the particle radius R P is the reference length, the velocity U of the approaching fluid relative to the particle is the reference velocity, the solubility of oil in the biofilm, c A,λ/β , and in the aqueous phase, c A,λ/υ , is the reference concentration for the respective phase. In particular, the following dimless quantities are defined The equation set defined in Equations (9a)-(9f) can be broken down into two subproblems that can be solved independently. The external mass transport problem defined by Equations (9d)-(9f) must be solved first, in order to determine the mass transfer coefficient k p/υ and the Biot number. As will be shown, the specific value of the solute concentration at the particle surface affects the concentration field in the υ-phase, but not the Biot number. Thereafter, the internal mass transport problem defined by Equations (9a)-(9c) must be solved in order to determine the overall dissolution rate at the surface of the oily core.

Advection-Dominated Transport in the Aqueous Phase without Bioreaction
In the absence of bioreaction (Da υ = 0), the external mass transport problem for the unbounded aqueous domain (Ω υ ) obtains the form and can be solved analytically in the limits of very low (Pe υ 1) or high Péclet number (Pe υ 1) [37,38]. Here, the high-Péclet regime is of primary interest and, thus, the derivation of the pertinent analytical solution is outlined. In spherical coordinates, for an axisymmetric concentration field (i.e., independent of the azimuthal angle), Equation (11a) obtains the detailed form v υ,r ∂c Aυ ∂r Moreover, for creeping Newtonian flow past a rigid sphere, the velocity components are [39] v υ,r (r, θ) = − 1 − 3 2r For advection-dominated mass transport, the change in the concentration from the value at the sphere surface (c Aυ = const.) to the bulk value away from the sphere (c Aυ = 0) is expected to occur within a thin boundary layer around the sphere. Upon this consideration, the following (dimless) independent variable is introduced to measure the distance from the sphere surface, within the boundary layer. On the basis that the thickness of the concentration boundary layer is small as compared to the radius of the sphere, i.e., y 1, the velocity terms can be simplified and certain diffusion terms can be neglected in Equation (12). In particular, order of magnitude analysis shows that the terms of tangential diffusion and normal diffusion due to surface curvature are much less important than the normal diffusion term. Under the boundary layer approximation, the final form of the reduced advection-diffusion equation is − 3 2 y 2 cos θ ∂c Aυ ∂y A detailed derivation of the above equation and the development of an analytical solution by means of a similarity transformation is given in [22] (pp. 80-87) and [39] (pp. 414-417). The exact solution of Equation (15) can be expressed as follows where C 2 is an integration constant given by and χ(y, θ) is a composite variable defined as with The concentration field given in Equation (16) is used to determine the diffusive mass flux and, ultimately, the average mass transfer rate from the particle surface to the aqueous phase The concentration derivative is calculated using the fundamental theorem of calculus, as follows Bioengineering 2018, 5, 15 9 of 29 and, after some operations, the final expression for the dissolution rate from the particle surface to the υ-phase, is given by the following expression where c Aυ ( R P ) = c A,λ/υ c Aβ (1) = H A,υ/β c Aβ ( R P ), and Here, k 0 p/υ is the mass transfer coefficient and the "0" superscript denotes the absence of bioreaction in the bulk aqueous phase. At this point, it is very useful to introduce the Sherwood number which is defined as follows and represents a dimless mass transfer coefficient. The above correlation underestimates the Sherwood number for about 10% for Pe υ > 100 [37] and, as expected, provides a wrong asymptotic value for Pe υ → 0 . By simply adding the value of the Sherwood number that corresponds to diffusion-only (i.e., Sh = 2 for Pe = 0), the following improved correlation is obtained Levich suggested the above superposition based on the rationale that the resistances to mass transfer by diffusion and advection act in parallel [22]. The estimates of Equation (22b) agree with numerical data within approximately 7% for the entire range of Pe. At this point, two remarks are in order. First, the mass transfer coefficient that appears in the Biot number does not depend on the, yet unknown, interfacial concentration of the solute. Second, the definition given in Equation (18) and the final expression given in Equation (20) for the surface averaged mass transfer rate were introduced earlier in the derivation of the modified boundary condition given in Equation (8).

Advection-Dominated Transport and Homogeneous Bioreaction in the Aqueous Phase
Following the analysis presented previously for advection dominated mass transport under the boundary layer theory approximation, the reduced form of the advection-diffusion-reaction equation with the same boundary conditions as given in Equations (11b)-(11c). To the best of our knowledge, an analytical solution is not available for the above partial differential equation. Approximate solutions have been developed using the empirical θ-expansion method of Yuge, perturbation analysis, and numerical methods [37,[40][41][42][43][44][45]. For engineering applications, the following simple correlation has been proposed for the Sherwood number [42,45] Sh p/υ ≡ k p/υ (2 R P ) with where Ha is the Hatta modulus and Sh 0 p/υ is the Sherwood number given by Equation (22) for the case of no bioreaction in the aqueous phase. The heuristic correlation given in Equation (24) is based on the film theory approximation and has been shown to provide an acceptable fit to more accurate numerical data. It is used here to provide estimates for the dissolution rate from the particle surface to the υ-phase, through the following relation

Diffusion and Reaction in the Biofilm Phase
The solution of the internal mass transfer problem is expressed as follows where h T is the Thiele modulus for homogeneous reaction in a spherical shell, and the integration constants are given by The concentration at the particle surface (υβ-interface) is now given by the expression and the concentration derivative at the core surface (λβ-interface) is The overall dissolution rate is where S λβ = 4π R 2 C is the area of the spherical core and the mass transfer coefficient is given by Again, it is convenient to define the Sherwood number The final quantity of interest is the volume averaged concentration of solute A in the β-phase which is later necessary in the determination of the particle size evolution. Here, V β = V P − V C is the volume of the biofilm shell, with V P = π D 3 P /6 and V C = π D 3 C /6. After some algebraic operations, the final expression for the volume averaged concentration is with where the concentration and its derivative are given in Equations (30) and (31), respectively.

Evolution of the Particle Size: Analysis of the Overall Mass Balances
The knowledge of the oil dissolution rate at the oil-biofilm and biofilm-water interfaces as well as of the average oil concentration in the biofilm can be used to determine the change in the dimensions of the compound particle over time. This is achieved through the analysis of the overall mass balance for the λand β-phases.

Overall Mass Balance for the λ-Phase
Upon considering the λ-phase as an open system that may exchange mass with the surrounding phases, the integral form of the mass balance is The term on the left hand side of the above equation represents the net accumulation of mass in the Ω λ -region. On the right hand side, the first term represents the change in the mass because of reaction in the Ω λ -region, and the second term represents the net influx of mass passing through the λβ-interface. Considering that the density of the λ-phase is constant, the accumulation term gives where D C,t = 2 R C,t is the diameter of the oily core at the t instant of time. The direct interfacial uptake of oil is modeled as a surface reaction occurring uniformly over the droplet surface and the reaction term obtains the form where δ λβ is Dirac's delta function concentrated on the λβ-interface, S λβ,t = π D 2 C,t is the surface area of the oily core, and r A,λβ is the constant reaction rate given in Equation (5).
The last term in Equation (38) represents the diffusive mass flux of oil across the λβ-interface. For completely immiscible phases, this term should be nil. For the problem at hand, the dissolution of the oil droplet is considered to be sufficiently slow ( k λ/β << U) so as not to have an appreciable impact on fluid dynamics but, nonetheless, this results in a non-zero diffusive flux across the droplet surface. Therefore, for this term we have with the overall dissolution rate given by Equation (32). Substitution of Equations (39)- (41) into where the Sherwood number varies along with the changing particle dimensions over time.

Overall Mass Balance for the β-Phase
The diameter of the compound particle also changes as the dissolving oily core is shrinking and its evolution is determined by the overall mass balance for the β-phase. We have that On the right hand side of the above equation, the second term could represent cell migration into the oily phase and the third term could represent the attachment of suspended cells to the biofilm or biofilm detachment and entrainment into the aqueous phase. However, for the problem at hand, both of these terms are considered to be nil. For the accumulation term on the left hand side, we have The rate of change in the biofilm mass, which is caused by the growth of cells and the synthesis of the extracellular matrix, is considered to be proportional to the microbial cell proliferation rate, i.e., Therefore, for the first term on the right hand side of Equation (43), we have where c Aβ t is the volume averaged concentration of oil in the β-phase at the t time instant, and is given by the expression in Equation (36). Substitution of Equations (44) and (45) into Equation (43), gives and, further, substitution of Equation (36) for the average concentration, after some operations, gives with

Compact and Dimless Forms of the Coupled ODEs
It is very convenient to express the coupled ordinary differential equations (ODEs) given in Equations (42) and (47), into the following compact form with Here, k srn is the droplet shrinking rate caused by direct interfacial uptake, k dis is the droplet shrinking rate caused by dissolution into the surrounding biofilm and aqueous phases, and k grt is the biofilm expansion rate due to growth.
The Damköhler and Thiele numbers that appear in the expressions given in Equations (34) and (37) for the Sherwood number Sh λ/β and the J C,t parameter, respectively, depend explicitly on the changing diameter of the compound particle. The situation might be a little bit more complex for the Péclet number in the case of a freely rising or sinking particle because the Stokes velocity also depends on the changing particle dimensions and density as follows where ∆ ρ t = | ρ υ − ρ P,t | is the excess density of the compound particle as compared to the density of the surrounding aqueous phase, ρ P,t = ϕ λ,t ρ λ + (1 − ϕ λ,t ) ρ β is the density of the compound particle, and ϕ λ,t = D 3 C,t / D 3 P,t is the volume fraction of the oily core. For a given set of parameters and initial conditions, the coupled ODEs given in Equations (49a) and (49b) can be solved numerically using, for instance, the explicit Euler or the classical Runge-Kutta method. In this work, both methods have been successfully implemented with an in-house Fortran code.
One step further, it is useful to establish a dimless form for the coupled ODEs that describe the evolution of the dimensions of the compound particle. For this purpose, a scaled characteristic diffusion time is introduced as follows Multiplication of both parts of Equations (49a) and (49b) with τ D / D P,0 gives with For the above dimless ODEs, it is only required to define dimless ratios (diffusivity, solubility, density) and the initial values of the dimless moduli. Thereafter, at each time instant the Damköhler and the Thiele moduli are updated as follows For the Péclet number, if the particle velocity is held constant we have Pe υ,t = D P,t Pe υ,0 , whereas if the particle is freely rising or sinking we have with the dimless excess density given by ∆ρ t = |1 − ρ P,t |/|1 − ρ P,0 |, and ρ P,t = ρ P,t / ρ υ . Two final remarks are in order. First, by setting δ β = 0, D C,t = D P,t , and c A,λ/β = c A,λ/υ , the compound particle model degenerates into a single-phase shrinking particle model. Second, the radius is the preferred characteristic length in the analysis of the local mass balances because it naturally arises with the spherical coordinate system. On the other hand, the particle diameter is used in the evolution of the particle dimensions because this length is determined experimentally in particle size analyses.

Results and Discussion
The main outcome of the theoretical model developed in Section 2 is the expression for the overall dissolution rate for a given particle configuration and the coupled ordinary differential equations for the evolution of the dimensions of the compound particle. In this section, the effect of key system parameters on the dissolution rate and the particle size evolution are investigated.

Overall Sherwood Number
According to Equation (32), for a given particle configuration and constant oil solubility, the overall dissolution rate increases with increasing Sherwood number. As already mentioned, the Sherwood number Sh λ/β represents the dimless mass transfer coefficient from the surface of the oily core to the surrounding biofilm shell and depends on the Biot number, the Thiele modulus, the thickness of the biofilm, and the solubility and diffusivity ratios. Among these parameters, the Biot number expresses the ratio of the external mass transfer rate (i.e., from the surface of the compound particle to the unbounded aqueous phase) over the characteristic intraparticle diffusion rate. The Biot number, in turn, depends on the Péclet and Damköhler numbers for the aqueous phase as well as on the solubility and diffusivity ratios. Figure 2 presents the dependence of the Biot number on the Péclet number for different values of the Damköhler number, keeping the other parameters constant. As expected, the Biot number increases with increasing Péclet and Damköhler numbers because the external mass transfer rate is enhanced by the contributions of advection and bulk bioreaction, respectively. For Pe υ = 0 and Da υ = 0, the solute moves away from the particle surface only by diffusion and the Biot number obtains the asymptotic value of Bi = H A,υ/β /Λ Aβ . For most solutes, diffusion within the biofilm is hindered by the extracellular matrix [46,47] and the diffusivity ratio is expected to be Λ Aβ ≤ 1. Furthermore, scarce experimental evidence suggests that the solubility of hydrophobic organic compounds might be significantly higher in the biofilm than in the aqueous phase (H A,υ/β ≤ 1). In the limit of exiguous solubility in the aqueous phase, i.e., H A,υ/β 1, the Biot number practically becomes nil and the dissolved oil is retained within the biofilm shell.
obtains the asymptotic value of Bi = , / Λ ⁄ . For most solutes, diffusion within the biofilm is hindered by the extracellular matrix [46,47] and the diffusivity ratio is expected to be Λ ≤ 1. Furthermore, scarce experimental evidence suggests that the solubility of hydrophobic organic compounds might be significantly higher in the biofilm than in the aqueous phase ( , / ≤ 1). In the limit of exiguous solubility in the aqueous phase, i.e., , / ≪ 1, the Biot number practically becomes nil and the dissolved oil is retained within the biofilm shell.    Figure 3a shows that the overall Sherwood number increases monotonically with increasing Biot and Thiele numbers, while keeping constant the other parameters. The Thiele number is a measure of the bioreaction rate over the diffusion rate within the biofilm shell. Faster bioreaction results in a steeper concentration gradient which, in turn, drives a higher rate of oil dissolution from the surface of the oily core. Of particular interest is the case of the vanishingly small Biot number which translates into the dissolved oil remaining trapped within the biofilm until complete biodegradation is achieved. Such a function would be of great practical importance and could perhaps be implemented by biofilms with hydrophobic or lipophilic biopolymers in their extracellular matrix [48]. This aspect deserves to be examined experimentally. biodegradation is achieved. Such a function would be of great practical importance and could perhaps be implemented by biofilms with hydrophobic or lipophilic biopolymers in their extracellular matrix [48]. This aspect deserves to be examined experimentally. Figure 3b presents a very interesting finding. If the Biot number is below a critical value (Bi ≈ 7 in this figure), then the overall Sherwood number increases monotonically as the dimless biofilm thickness increases. On the other hand, if the Biot number is above the critical value, then the Sherwood number decreases as the biofilm thickness increases from zero up to a critical value , and, beyond that value, the Sherwood number re-increases with increasing biofilm thickness. In order to elucidate this behavior, a detailed examination of the expression given in Equation (34) for the overall Sherwood number is required. It has been established that The first contribution contains the intertwined effects of transport and bioreaction in the biofilm and aqueous phases. With increasing biofilm thickness, this term increases or decreases monotonically up or down to the asymptotic value of 2Λ ℎ , depending on the relative importance of the internal (biofilm) and external (aqueous) resistances to mass transport. If the internal resistance to mass transport is lower than the external resistance (i.e., sufficiently low Biot and high Thiele numbers), then the first contribution increases up to the asymptotic value as the biofilm thickness increases. In the opposite case of higher internal resistance (i.e., sufficiently high Biot and low Thiele numbers), the first contribution decreases down to the asymptotic value with increasing biofilm thickness.
The second term in the formula for the overall Sherwood represents the effect of curvature on the concentration gradient that is evaluated at the core surface and increases monotonically with increasing biofilm thickness. In particular, each point on the oily core surface is projected to a surface element of finite area on the outer surface of the compound particle. As the distance between the two concentric spherical surfaces increases, the degree of geometric expansion also increases and causes a dilution in the solute concentration at the outer surface which, in turn, results in a higher concentration gradient. This geometric effect does not exist for a flat surface configuration. The interplay between the two contributions produces the pattern shown in Figure 3b.     figure), then the overall Sherwood number increases monotonically as the dimless biofilm thickness increases. On the other hand, if the Biot number is above the critical value, then the Sherwood number decreases as the biofilm thickness increases from zero up to a critical value δ β,crit and, beyond that value, the Sherwood number re-increases with increasing biofilm thickness. In order to elucidate this behavior, a detailed examination of the expression given in Equation (34) for the overall Sherwood number is required. It has been established that The first contribution contains the intertwined effects of transport and bioreaction in the biofilm and aqueous phases. With increasing biofilm thickness, this term increases or decreases monotonically up or down to the asymptotic value of 2Λ Aβ h T , depending on the relative importance of the internal (biofilm) and external (aqueous) resistances to mass transport. If the internal resistance to mass transport is lower than the external resistance (i.e., sufficiently low Biot and high Thiele numbers), then the first contribution increases up to the asymptotic value as the biofilm thickness increases. In the opposite case of higher internal resistance (i.e., sufficiently high Biot and low Thiele numbers), the first contribution decreases down to the asymptotic value with increasing biofilm thickness.
The second term in the formula for the overall Sherwood represents the effect of curvature on the concentration gradient that is evaluated at the core surface and increases monotonically with increasing biofilm thickness. In particular, each point on the oily core surface is projected to a surface element of finite area on the outer surface of the compound particle. As the distance between the two concentric spherical surfaces increases, the degree of geometric expansion also increases and causes a dilution in the solute concentration at the outer surface which, in turn, results in a higher concentration gradient. This geometric effect does not exist for a flat surface configuration. The interplay between the two contributions produces the pattern shown in Figure 3b.
The critical Biot number, above which the biofilm shell acts as a diffusive barrier and hinders the transport of oil compounds from the surface of the core to the surrounding aqueous phase, can be determined by the following condition With reference to Figure 3b, the above condition states that the constant Biot curve for the critical Biot value is normal to the Sh λ/β axis with abscissa δ β = 0. The first derivative of the Sherwood number with respect to the dimless biofilm thickness is determined from Equation (34) and, after some operations, obtains the following form Substitution of the above expression in Equation (57), gives the following expression for the critical Biot The critical biofilm thickness, below which the biofilm hinders mass transport, corresponds to the abscissa of the minimum in the constant Biot curves (for Bi > Bi crit ). Therefore, for given Biot and Thiele numbers, it can be determined as the root of the nonlinear algebraic equation that is obtained by setting the first derivative given in Equation (58) equal to zero. Another straightforward way to obtain an estimate for the critical biofilm thickness is to consider that at this value the first contribution in the Sherwood expression is sufficiently close to the asymptotic value of 2Λ Aβ h T . Therefore, by demanding that tanh h T δ β,crit = 0.99, the following simple estimate is obtained. The critical biofilm thickness depends also on the Biot number but, as can be seen in Figure 3b, the dependence is weak and, thus, the above estimate is sufficient for practical purposes.

Relative Importance of the Bioreaction and Dissolution Processes
Part of the oil that dissolves at the oil-biofilm interface is biodegraded within the biofilm and the rest is released into the water column; where it might, or might not, be further degraded by suspended microbes. A key issue concerns the bioreactive effectiveness of the biofilm shell. The amount of dissolved oil that ends up either biodegraded or released can be determined by the overall mass balances (Section 2.4).
The overall dissolution rate W A,λ/β of oil at the surface of the oily core is given in Equation (32). The oil dissolution rate from the particle surface to the surrounding aqueous phase is given in Equation (26) and can be expressed in the following equivalent form using the relation c Aυ R P = c Aβ (1) c A,λ/υ = c Aβ (1)H A,υ/β c A,λ/β , with the dimless concentration c Aβ (1) given in Equation (30). The rate of oil bioreaction within the biofilm can be expressed as with the average concentration of oil in the biofilm shell given by Equation (36). The mass fractions of biodegraded oil in the biofilm and released oil in the water column are defined as respectively, with Φ brn + Φ dis = 1. Figure 4 presents the effect of the Thiele modulus on the biodegraded and released oil fractions for different values of the Biot number and the thickness of the biofilm shell. It is observed that the mass fraction of biodegraded oil increases with increasing Thiele modulus, decreasing Biot number, and increasing biofilm thickness. As a consequence, biofilms composed of fast oil-degrading microbes and lipophilic extracellular matrix would be ideal for retaining and biodegrading oil compounds in practical applications.
respectively, with Φ + Φ = 1. Figure 4 presents the effect of the Thiele modulus on the biodegraded and released oil fractions for different values of the Biot number and the thickness of the biofilm shell. It is observed that the mass fraction of biodegraded oil increases with increasing Thiele modulus, decreasing Biot number, and increasing biofilm thickness. As a consequence, biofilms composed of fast oil-degrading microbes and lipophilic extracellular matrix would be ideal for retaining and biodegrading oil compounds in practical applications.

Impact of the Péclet and Thiele Numbers on the Particle Size Evolution
The evolution of the dimensions of the compound particle is determined by all the factors that affect the dissolution of oil into the surrounding phases, the direct uptake of oil at the surface of the oily core, and the volumetric growth of the biofilm phase. First, we examine the effects of the bioreaction in the biofilm (expressed by the Thiele number) and of the particle velocity (expressed by the Péclet number), while considering that the rates of direct uptake and biofilm growth are nil. It is very convenient to use the dimless form of the coupled ODEs given in Equations (53a) and (53b), as it is only required to define certain dimless quantities without specifying the values of solubilities, kinetic and other system parameters. Figure 5 presents the strong effect of the initial Thiele number on the evolution of the particle dimensions, while keeping all other parameters constant. As expected, higher Thiele numbers result in higher shrinking rates and faster consumption of the oily core. An interesting feature is that the temporal change in the dimensions of the particle is non-linear. For a given Thiele number, the diameter of the oily core decreases with an increasing shrinking rate (Figure 5a, concave function), whereas the diameter of the compound particle decreases with a decreasing shrinking rate (Figure 5b, convex function) until reaching an asymptotic value that corresponds to a residual particle that contains only biofilm.

Impact of the Péclet and Thiele Numbers on the Particle Size Evolution
The evolution of the dimensions of the compound particle is determined by all the factors that affect the dissolution of oil into the surrounding phases, the direct uptake of oil at the surface of the oily core, and the volumetric growth of the biofilm phase. First, we examine the effects of the bioreaction in the biofilm (expressed by the Thiele number) and of the particle velocity (expressed by the Péclet number), while considering that the rates of direct uptake and biofilm growth are nil. It is very convenient to use the dimless form of the coupled ODEs given in Equations (53a) and (53b), as it is only required to define certain dimless quantities without specifying the values of solubilities, kinetic and other system parameters. Figure 5 presents the strong effect of the initial Thiele number on the evolution of the particle dimensions, while keeping all other parameters constant. As expected, higher Thiele numbers result in higher shrinking rates and faster consumption of the oily core. An interesting feature is that the temporal change in the dimensions of the particle is non-linear. For a given Thiele number, the diameter of the oily core decreases with an increasing shrinking rate (Figure 5a, concave function), whereas the diameter of the compound particle decreases with a decreasing shrinking rate (Figure 5b, convex function) until reaching an asymptotic value that corresponds to a residual particle that contains only biofilm.  Figure 6 presents the effect of the initial Péclet number on the evolution of the particle dimensions, for the case of an inert biofilm shell (non-reactive, non-growing) and while keeping all other parameters constant. Upon the assumption that the diffusion coefficient and the initial particle size are held constant, the different Péclet numbers are associated with different characteristic velocities. Two cases are considered for the mode of change in the characteristic velocity as the particle shrinks. In the "free rise/fall" mode, the particle is considered to rise or fall in the water column under the action of gravity with the velocity given by Stokes' formula in Equation (51) and the Péclet number updated according to Equation (56). In the "const. U" mode, the particle is considered to be carried by the aqueous stream at constant velocity. First of all, as expected, a higher Péclet number leads to faster dissolution and shrinkage of the oily core. Furthermore, for a given initial Péclet number, the velocity mode appears to have an appreciable effect on the shrinking rate. In particular, the "free rise/fall" mode exhibits an interesting dependence on the biofilm density ( Figure 7). Typically, the density of the oily phase is lower than that of the surrounding aqueous phase. Therefore, oil droplets without a biofilm shell tend to rise when released in a water column. The situation is different for compound droplets because the biofilm density shows significant variability as it depends on the content and type of cells and extracellular polymers. If the biofilm density is lower than or similar to the density of the aqueous phase, then the compound particle rises with decreasing velocity and Péclet number as the oily core is consumed over time (cases of ≤ 1 in Figure 7). On the other hand, if the biofilm density is higher than the density of the aqueous phase, then the compound particle rises until it becomes neutrally buoyant and, thereafter, begins to sink with increasing velocity as the oily core is shrinking (cases of > 1 in Figure 7).  Figure 6 presents the effect of the initial Péclet number on the evolution of the particle dimensions, for the case of an inert biofilm shell (non-reactive, non-growing) and while keeping all other parameters constant. Upon the assumption that the diffusion coefficient and the initial particle size are held constant, the different Péclet numbers are associated with different characteristic velocities. Two cases are considered for the mode of change in the characteristic velocity as the particle shrinks. In the "free rise/fall" mode, the particle is considered to rise or fall in the water column under the action of gravity with the velocity given by Stokes' formula in Equation (51) and the Péclet number updated according to Equation (56). In the "const. U" mode, the particle is considered to be carried by the aqueous stream at constant velocity. First of all, as expected, a higher Péclet number leads to faster dissolution and shrinkage of the oily core. Furthermore, for a given initial Péclet number, the velocity mode appears to have an appreciable effect on the shrinking rate. In particular, the "free rise/fall" mode exhibits an interesting dependence on the biofilm density ( Figure 7). Typically, the density of the oily phase is lower than that of the surrounding aqueous phase. Therefore, oil droplets without a biofilm shell tend to rise when released in a water column. The situation is different for compound droplets because the biofilm density shows significant variability as it depends on the content and type of cells and extracellular polymers. If the biofilm density is lower than or similar to the density of the aqueous phase, then the compound particle rises with decreasing velocity and Péclet number as the oily core is consumed over time (cases of ρ β ≤ 1 in Figure 7). On the other hand, if the biofilm density is higher than the density of the aqueous phase, then the compound particle rises until it becomes neutrally buoyant and, thereafter, begins to sink with increasing velocity as the oily core is shrinking (cases of ρ β > 1 in Figure 7).

Impact of Biofilm Growth and Direct Uptake on the Particle Size Evolution
The effect of a growing biofilm on the temporal evolution of the particle configuration is somewhat intricate as a thicker biofilm might not always enhance the droplet shrinking and oil biodegradation rates. If , / ⁄ < 1.2, that is if the oil is more soluble and mobile in the biofilm than in the aqueous phase, then the internal resistance to mass transport is also lower than the external resistance (Bi < Bi ) and a net increase in the amount of biofilm due to growth results in higher rates of oil dissolution and droplet shrinking (Figure 8). Under such conditions, the enhancement of the biodegradation process is more profound with thicker biofilms, i.e., for higher values of the biofilm yield coefficient / . In Figures 8 and 9, the value of / = 0 corresponds to the case where the initial amount of biofilm remains constant over time and is just redistributed around the oily core as the latter shrinks. For / > 0, additional volume of biofilm is produced.

Impact of Biofilm Growth and Direct Uptake on the Particle Size Evolution
The effect of a growing biofilm on the temporal evolution of the particle configuration is somewhat intricate as a thicker biofilm might not always enhance the droplet shrinking and oil biodegradation rates. If , / ⁄ < 1.2, that is if the oil is more soluble and mobile in the biofilm than in the aqueous phase, then the internal resistance to mass transport is also lower than the external resistance (Bi < Bi ) and a net increase in the amount of biofilm due to growth results in higher rates of oil dissolution and droplet shrinking (Figure 8). Under such conditions, the enhancement of the biodegradation process is more profound with thicker biofilms, i.e., for higher values of the biofilm yield coefficient / . In Figures 8 and 9, the value of / = 0 corresponds to the case where the initial amount of biofilm remains constant over time and is just redistributed around the oily core as the latter shrinks. For / > 0, additional volume of biofilm is produced. Figure 7. Impact of the biofilm density on the evolution of the Péclet number for a shrinking compound particle that freely rises/falls in a water column. The compound particle becomes neutrally buoyant at the time instants indicated by the arrows. The parameter values are: Pe υ,0 = 100; δ β,0 = 0.1; h T = 0;

Impact of Biofilm Growth and Direct Uptake on the Particle Size Evolution
The effect of a growing biofilm on the temporal evolution of the particle configuration is somewhat intricate as a thicker biofilm might not always enhance the droplet shrinking and oil biodegradation rates. If H A,υ/β /Λ Aβ < 1.2, that is if the oil is more soluble and mobile in the biofilm than in the aqueous phase, then the internal resistance to mass transport is also lower than the external resistance (Bi < Bi crit ) and a net increase in the amount of biofilm due to growth results in higher rates of oil dissolution and droplet shrinking (Figure 8). Under such conditions, the enhancement of the biodegradation process is more profound with thicker biofilms, i.e., for higher values of the biofilm yield coefficient Y β/A . In Figures 8 and 9, the value of Y β/A = 0 corresponds to the case where the initial amount of biofilm remains constant over time and is just redistributed around the oily core as the latter shrinks. For Y β/A > 0, additional volume of biofilm is produced. On the other hand, if , / ⁄ > 1.2, then the internal resistance to mass transport is higher than the external resistance (Bi > Bi ) and a net increase in the amount of biofilm due to growth results in lower rates of oil dissolution and droplet shrinking as compared to the non-growing case ( Figure 9). The peculiar trend observed in Figure 9b for the dissolution rate is attributed to the role of biofilm as a diffusive barrier, which is discussed in detail in Section 3.1. In short, the dissolution rate is defined in Equation (54) as ( ) = 2 ℎ / , ⁄ and follows closely the dependence of the Sherwood number on the dimless biofilm thickness. As the oily core shrinks, the dimless biofilm thickness increases and, for Bi > Bi , the Sherwood number decreases. This occurs until the biofilm exceeds the critical biofilm thickness, > , , so as the enhancement caused by the curvature effect to supersede the attenuation caused by the diffusive barrier. Thereafter, the Sherwood number increases with increasing biofilm thickness (see also Figure 3b and the discussion in Section 3.1).
(a) (b)  On the other hand, if , / ⁄ > 1.2, then the internal resistance to mass transport is higher than the external resistance (Bi > Bi ) and a net increase in the amount of biofilm due to growth results in lower rates of oil dissolution and droplet shrinking as compared to the non-growing case ( Figure 9). The peculiar trend observed in Figure 9b for the dissolution rate is attributed to the role of biofilm as a diffusive barrier, which is discussed in detail in Section 3.1. In short, the dissolution rate is defined in Equation (54) as ( ) = 2 ℎ / , ⁄ and follows closely the dependence of the Sherwood number on the dimless biofilm thickness. As the oily core shrinks, the dimless biofilm thickness increases and, for Bi > Bi , the Sherwood number decreases. This occurs until the biofilm exceeds the critical biofilm thickness, > , , so as the enhancement caused by the curvature effect to supersede the attenuation caused by the diffusive barrier. Thereafter, the Sherwood number increases with increasing biofilm thickness (see also Figure 3b and the discussion in Section 3.1).
(a) (b)  On the other hand, if H A,υ/β /Λ Aβ > 1.2, then the internal resistance to mass transport is higher than the external resistance (Bi > Bi crit ) and a net increase in the amount of biofilm due to growth results in lower rates of oil dissolution and droplet shrinking as compared to the non-growing case ( Figure 9). The peculiar trend observed in Figure 9b for the dissolution rate is attributed to the role of biofilm as a diffusive barrier, which is discussed in detail in Section 3.1. In short, the dissolution rate is defined in Equation (54) as k dis (τ) = 2Sh λ/β /D P,t and follows closely the dependence of the Sherwood number on the dimless biofilm thickness. As the oily core shrinks, the dimless biofilm thickness increases and, for Bi > Bi crit , the Sherwood number decreases. This occurs until the biofilm exceeds the critical biofilm thickness, δ β > δ β,crit , so as the enhancement caused by the curvature effect to supersede the attenuation caused by the diffusive barrier. Thereafter, the Sherwood number increases with increasing biofilm thickness (see also Figure 3b and the discussion in Section 3.1). Figure 10 illustrates the potential effect of the biodegradation and dissolution processes on the shrinking rate of the oily core. The parameter values are typical for the applications under consideration (see a detailed discussion in the next section), and have also been selected so as to exemplify that each mechanism might have a significant impact on the overall process. In real world applications, one or two or all mechanisms might act in parallel. For the direct uptake mechanism, the value of k srn = 10 is obtained in Equation (54) by setting the direct uptake rate k srn = 0.2 nm/s (Table 1), the initial droplet diameter D P,0 = 100 µm, the oil density ρ λ = 0.85 g/cm 3 , the oil diffusivity in water D Aυ = 10 −6 cm 2 /s, and the oil solubility in biofilm c A,λ/β = 17 µg/cm 3 .
Bioengineering 2018, 5, x FOR PEER REVIEW 21 of 28 Figure 10 illustrates the potential effect of the biodegradation and dissolution processes on the shrinking rate of the oily core. The parameter values are typical for the applications under consideration (see a detailed discussion in the next section), and have also been selected so as to exemplify that each mechanism might have a significant impact on the overall process. In real world applications, one or two or all mechanisms might act in parallel. For the direct uptake mechanism, the value of = 10 is obtained in Equation (54) by setting the direct uptake rate = 0.2 nm s ⁄ (Table 1), the initial droplet diameter , = 100 μm, the oil density = 0.85 g cm ⁄ , the oil diffusivity in water = 10 cm s ⁄ , and the oil solubility in biofilm ̃ , / = 17 μg cm ⁄ .

Implications for the Biodegradation of Crude Oil Microdroplets in the Sea
At this point, a naturally arising question concerns the values of the characteristic dimless moduli and other system parameters for real world applications. First of all, as mentioned in the introduction, the microdroplet might be rising, sinking or drifting along underwater sea currents. For light crude oil microdroplets with a diameter in the range of , = 10 − 100 μm and a density of = 0.85 g cm ⁄ , freely rising due to buoyancy through an aqueous water column with density = 1.02 g cm ⁄ and dynamic viscosity = 0.01 g (cm • s) ⁄ , the Stokes velocity given by Equation (51) is in the range of , ≅ 9 • 10 − 0.09 cm s ⁄ and the corresponding radius-based Reynolds number is Re ≈ 5 • 10 − 5 • 10 . The density of the biofilm is expected to be similar to or larger than the density of the aqueous phase, depending on the type and volume fraction of cells and extracellular biopolymers within the biofilm. For instance, the density of marine snow particles

Implications for the Biodegradation of Crude Oil Microdroplets in the Sea
At this point, a naturally arising question concerns the values of the characteristic dimless moduli and other system parameters for real world applications. First of all, as mentioned in the introduction, the microdroplet might be rising, sinking or drifting along underwater sea currents. For light crude oil microdroplets with a diameter in the range of D P,0 = 10 − 100 µm and a density of ρ λ = 0.85 g/cm 3 , freely rising due to buoyancy through an aqueous water column with density ρ υ = 1.02 g/cm 3 and dynamic viscosity µ υ = 0.01 g/(cm·s), the Stokes velocity given by Equation (51) is in the range of U S,0 ∼ = 9·10 −4 − 0.09 cm/s and the corresponding radius-based Reynolds number is Re υ ≈ 5·10 −5 − 5·10 −2 . The density of the biofilm is expected to be similar to or larger than the density of the aqueous phase, depending on the type and volume fraction of cells and extracellular biopolymers within the biofilm. For instance, the density of marine snow particles collected from the Gulf of Mexico after the Deepwater Horizon incident exhibited great variability with values of excess density ranging from 0.07 g/cm 3 to 0.36 g/cm 3 [49]. Therefore, for compound particles with a diameter in the range of 10 − 100 µm and density ρ P = 1.25 g/cm 3 , the settling Stokes velocity is also on the order of 10 −3 − 10 −1 cm/s. In the case of (almost) neutrally buoyant particles, the characteristic velocity is determined by the velocity of the underwater sea current that carries the particles. The velocity of sea currents varies over several orders with a magnitude of a few mm/s for the vertical velocity [50], several cm/s for the horizontal velocity in deep sea (depth larger than 300 m) [6,[50][51][52], and from tenths of cm/s up to a few m/s for the horizontal velocity near the sea surface [50][51][52]. For example, in the Gulf of Mexico, values of 2 − 3 mm/s have been measured for the vertical velocity [50] and an average value of 7.8 cm/s has been reported for the horizontal velocity at a depth of 1100 m [6]. Thus, for compound particles drifting along an underwater current with a velocity in the range of 0.1 − 10 cm/s, the Reynolds number is Re υ ≈ 5·10 −3 − 5.
At atmospheric conditions, the interfacial tension between crude oil and water is in the range of γ λυ = 10 − 30 dyn/cm, depending on the detailed composition of the two phases [53,54]. An even higher tension might be expected between a compound particle and water, especially if the biofilm that covers the oily core is strongly hydrophobic. For instance, it has been recently found that microbes of the species Bacillus subtilis secrete a hydrophobic protein called BslA, which accumulates in the outer layer of the biofilm and results in strong repellence of aqueous drops (contact angle > 90 • ) [48,55]. Nonetheless, to the best of our knowledge, specific values of the interfacial tension for biofilm-water systems have not been reported yet. For a characteristic velocity in the range of 10 −3 − 10 cm/s and an interfacial tension of 20 dyn/cm the capillary number is Ca = 5·10 −7 − 5·10 −3 . Furthermore, for compound particles with a diameter in the range of 10 − 100 µm and density of 1.25 g/cm 3 the Bond number is Bo = 3·10 −6 − 3·10 −4 . Finally, for a characteristic velocity in the range of 10 −3 − 10 cm/s, and for dissolved oil components with a diffusion coefficient on the order of D Aυ = 10 −6 cm 2 /s, the Péclet number is Pe υ = 0.5 − 5·10 4 . In view of the above data, it is concluded that the hypotheses set out in Section 2.2 for the hydrodynamics and solute transport problems are reasonable for most of the cases considered in this work. Future extension of the CPM formulation to low-but-finite Reynolds numbers and retention of all terms in the solute mass balance is expected to improve substantially the domain of validity of the proposed model. With regard to the bioreaction kinetic parameters, a major issue is raised. A theoretical model, focused on a specific spatial scale, requires the input of reaction rates and parameters which precisely correspond to the scale of focus. Typically, in oil biodegradation experiments, the physical state of the oil is not taken into explicit account and an apparent biodegradation rate is determined in terms of an average oil concentration that lumps together all forms of oil, i.e., dissolved, micellar, and/or dispersed in droplets. This lumped concentration and its spatial and temporal derivatives differ from the concentration field that is detected by the microbes in their microenvironment [56]; in the present context, from the concentrated oil detected by flatlanders at the oil-water interface, the concentration c Aυ detected by drifters in the bulk aqueous phase, and the concentration c Aβ detected by biofilm formers within the biofilm. Here, the term "apparent" is used to denote a reaction rate that pertains to a representative elementary volume with dimensions much larger than the size of a single microdroplet or a single microbial cell. For systems with microscale heterogeneity, such as porous media and multiphase dispersions, the apparent reaction rate incorporates the effects of the microscale structure and transport mechanisms and is, therefore, lower than the intrinsic (transport-free) reaction rate [57]. This important issue has recently attracted the attention of researchers working on the biodegradation of crude oil [58][59][60][61][62][63][64]. For example, it has been nicely demonstrated that the apparent biodegradation rate increases with decreasing droplet size, while keeping all other system parameters constant [30,58,61]. However, existing kinetic data for apparent reaction rates of oily compounds are not consistent with the CPM formulation, and can only be used to obtain lower bounds for the Damköhler and Thiele numbers.
Recent studies on the biodegradation of crude oil that was released by the Macondo well in the Gulf of Mexico after the Deepwater Horizon blowout, report that the apparent half-life of many biodegradable hydrocarbons is in the range of T 1/2 = 0.1 − 10 d (at~5 • C) [7,30,62,63]. The first-order reaction constant is related to the half-life as k 1,α = ln 2/ T 1/2 and obtains values in the range of k 1,α ≈ 0.7 − 7d −1 . Therefore, for microdroplets of diameter D P,0 = 100 µm and water diffusivity D Aυ = 10 −6 cm 2 /s, the Damköhler number for the aqueous phase obtains values in the range of Da υ = 2·10 −5 − 2·10 −3 . The diffusion coefficient of oil compounds in the interstitial space of biofilms is expected to be lower, by one or more orders of magnitude, than the water diffusivity [46,47] and, thus, the Thiele number for the biofilm phase is in the range of h T = 0.0045 − 0.045 (assuming Λ Aβ = 0.1). According to the CPM formulation, Equation (4), the first-order reaction constant depends not only on the kinetic parameters ( µ m,α , K S,α , Y C/A,α ), but also on the concentration B α of active microbial cells. In natural ecosystems, the concentration of cells residing in biofilms might be several orders higher than the concentration of suspended cells. For example, after the Deepwater Horizon incident, the concentration of marine microbes in the underwater oil plume was B υ = 10 4 − 10 6 cells/cm 3 [7], whereas the cell concentration might reach values of B β = 10 8 − 10 10 cells/cm 3 within marine snow and biofilms, depending on the size of individual cells and the cell volume fraction. Vilcáez et al. [32] developed a theoretical model for the biodegradation of droplet populations on the basis of a shrinking particle model that accounts only for the direct uptake mode. They also analyzed kinetic data from the literature and reported lumped parameters for three groups of hydrocarbons, namely alkanes, monoaromatics (BTEX), and polyaromatics (PAHs). For 100 µm-sized droplets, the data of Vilcáez et al. (Table 1) give values in the range of Da υ = 0.0001 − 0.0003 for the Damköhler number in the aqueous phase and h T = 1.02 − 1.96 for the Thiele number in the biofilm phase. Consequently, based on the available data, bioreaction is expected to be of considerable importance within the biofilm phase, whereas it is dominated by advection and diffusion in the aqueous phase. Nonetheless, because of the significant uncertainty with regard to the consistency between the CPM formulation and existing kinetic data for oily substrates, this discussion must be extended once data from microscale experiments become available.

Conclusions
In this paper, a compound particle model of the core-shell type is developed for the microbial degradation of solitary oil microdroplets and takes into account three fundamental biodegradation modes, namely the direct interfacial uptake at the oil surface, the bioreaction in the bulk aqueous phase, and the bioreaction in a biofilm formed around the droplet. Previous relevant models account only for the direct uptake mode. The major results of the theoretical analysis include an expression for the overall dissolution rate for a given particle configuration and two coupled ordinary differential equations for the evolution of the dimensions of the compound particle. An interesting finding is that biofilms consisting of a high concentration of fast oil-degrading microbes and lipophilic biopolymers (corresponding to a low-Biot and high-Thiele regime) are expected to be ideal for oil biodegradation applications because they retain the dissolved oil until complete degradation, instead of releasing it into the water column. The model is based on a large set of simplifying, yet justifiable, hypotheses; most of which rely on the consideration of microsized droplets with an immobilized interface by the presence of microbes and biopolymers. One of the most important hypotheses in the model is that the compound particle moves like a rigid sphere. This hypothesis negates the need to define the biofilm mechanics, which might range from a fluid-like to a solid-like behavior, depending on the composition of the biofilm and the applied stresses [65][66][67][68]. Another very important hypothesis is that the oily phase is treated as a single component, whereas crude oil and most natural or artificial oils are multi-component mixtures. Selective or faster biodegradation of certain components (e.g., alkanes) will result in concentration gradients within the oil droplet. The effects of these gradients on the droplet shrinking rate is expected to be minimal only if intra-droplet diffusion is much faster than diffusion in the surrounding biofilm phase (i.e., D Aλ / D Aβ 1). Besides the development of concentration gradients, selective biodegradation will also change the mass fraction of each oil compound within the droplet. In turn, the density of the droplet, ρ λ , will also vary with time. For instance, if light alkanes are consumed faster than heavier compounds like PAHs, then the velocity of a compound particle undergoing free rise will decrease faster and the impact on the temporal evolution of the Péclet number might be appreciably stronger than that shown in Figure 7. These hypotheses will be relaxed in future work by adding more physical and mathematical complexity into the model formulation.