Numerical study of plume patterns in the chemotaxis–diffusion–convection coupling system

Numerical study of plume patterns in the chemotaxis–diffusion–convection coupling system

Yannick Deleuze Chen-Yu Chiang Department of Engineering Science and Ocean Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei, Taiwan Marc Thiriet Tony W.H. Sheu

A chemotaxis–diffusion–convection coupling system for describing a form of buoyant convection in which the fluid develops convection cells and plume patterns will be investigated numerically in this study. Based on the two-dimensional convective chemotaxis-fluid model proposed in the literature, we developed an upwind finite element method to investigate the pattern formation and the hydrodynamical stability of the system. The numerical simulations illustrate different predicted physical regimes in the system. In the convective regime, the predicted plumes resemble Bénard instabilities. Our numerical results show how structured layers of bacteria are formed before bacterium rich plumes fall in the fluid. The plumes have a well defined spectrum of wavelengths and have an exponential growth rate, yet their position can only be predicted in very simple examples. In the chemotactic and diffusive regimes, the effects of chemotaxis are investigated. Our results indicate that the chemotaxis can stabilize the overall system. A time scale analysis has been performed to demonstrate that the critical taxis Rayleigh number for which instabilities set in depends on the chemotaxis head and sensitivity. In addition, the comparison of the differential systems of chemotaxis–diffusion–convection, double diffusive convection, and Rayleigh-Bénard convection establishes a set of evidences that even if the physical mechanisms are different at the same time the dimensionless systems are strongly related to each other.

Keywords: Chemotaxis, Finite element method, Incompressible viscous fluid flow, Coupled KS–NS

concentration of oxygen Greek symbols
concentration of oxygen in air solute expansion coefficient
diffusivity thermal expansion coefficient
growth rate of the plume amplitude bacterial oxygen consumption rate
container height dimensionless container width
chemotaxis head dynamic viscosity
vertical unit vector upwards kinematic viscosity
Le Lewis number fluid density
areal number density of bacteria volumetric mass density of a bacterium
unit outward normal vector
initial number density of bacteria Superscripts
initial average number density of bacteria dimensionless quantities
Pr Prandtl number Subscripts
Ra Rayleigh number bacterium
dimensionless chemotaxis sensitivity critical
dimensional chemotaxis sensitivity convection
solute concentration diffusion
time mass
time scale of bacterium transport oxygen
temperature solute
velocity vector taxis
velocity in x-direction
velocity in y-direction
volume of a bacterium
coordinate axes
Table 1: Nomenclature description

1 Introduction

Taxis refers to the collective motion of cells or an organism in response to an attractant gradient. The nature of the attractant stimulus can be of chemical (chemotaxis), physical (baro-, electro-, magneto-, phono-, photo-, and thermotaxis), or mechanical (hapto- and rheotaxis) origins. Chemotaxis refers to cell movement primed by a external chemical signal that can either be emitted by the same population of cells or created by an external source [1]. In particular, aerotaxis is related to the movement toward a gradient of increasing oxygen concentration.

The phenomenon that couples chemotaxis, diffusion, and convection has been illustrated by experiments on suspensions of bacteria in a container filled with water [2, 3]. Oxygen diffuses in the container from the surface. As bacteria consume oxygen, oxygen concentration falls everywhere except at the surface, hence creating a vertical concentration difference. Bacteria move up to higher concentration of oxygen and quickly get densely packed below the surface in a relatively thin layer. Subsequently, Rayleigh–Bénard-like instability appears near the surface. These dynamical instabilities exhibit complex convection patterns with plumes of bacteria falling in the fluid.

Chemotaxis–diffusion–convection is a particular case of the so-called bioconvection. Bioconvection is a more general term for suspensions of swimming microorganisms which are denser than the solvent fluid. Bioconvection and the different mechanisms of upswimming have been reviewed [4].

Mathematical modeling of chemotaxis has been introduced by authors in [5] and [6]. Mathematical analysis mainly focused on pattern formation of microorganisms and blow-up phenomena in finite time. Key contributions on chemotactic collapse have been brought by authors in [7, 8, 9, 10, 11, 12]. Similarly, angiogenesis is the motion of endothelial cells to form new blood vessels from preexisting vessels to locally supply oxygen. During tumor growth, tumor cells secrete a set of substances to attract endothelial cells [13, 14]. In angiogenesis models, chemotactic collapse was not observed [15, 16]. In this study, the model equations (1-4) are used to describe chemotactic response of bacterial suspensions. The unstable agglomeration of bacterial cells at the surface leads to a descent of cell-rich plume patterns together with the formation of high-speed jets between counter-rotating vortices.

Formation and stability of plumes result from the balance between chemotaxis, diffusion, and convection of bacteria. The particular impact of each mechanism still needs to be understood. Chemotaxis is known to bring instability in system and leads to aggregation, but may also have a stabilizing effect.

The linear stability analysis of the chemotaxis-diffusion-convection system showed that a condition for linear instability depends on the taxis Rayleigh number [2]. The taxis Rayleigh number is defined as the ratio of buoyancy and viscosity forces times the ratio of momentum and cell diffusivity (table 2). Below a critical value (), then the process remains stable. From experiments, several stages were observed starting from the upward bacteria accumulation and leading to hydrodynamic formation of plumes. A weakly nonlinear stability analysis was conducted to investigate the stability of hexagon and roll patterns formed by the system of equations (1-4) [17]. The hydrodynamic vortices formed by convection strengthen circulation of fluid and enhance the intake of oxygen into the solvent [18]. Global existence for the chemotaxis–Stokes system for small initial bacterial population density was proved in [19]. Then, global existence for the chemotaxis–Navier-Stokes system for a large initial bacterial population density as well as global existence of 3D weak solutions for the chemotaxis-Stokes equations were proved in [20].

Table 2: Representative dimensionless numbers involved in the double diffusive convection (DDC), chemotaxis–diffusion–convection (CDC), and Rayleigh-Bénard convection (RBC). Subscripts , , , , , stand for mass, thermal, taxis, solute, oxygen, and bacterium, respectively; is the acceleration due to gravity; and are the coefficients for thermal and solute expansion, respectively; is the kinematic viscosity.

A detailed numerical study of the plume formation and merging that was related to the convergence of Rayleigh-Bénard-type patterns was carried out in [21]. The shape and number of plumes can be controlled by initial bacterial population density. However, the sites of plumes have not been predicted. The convergence toward numerically stable stationary plumes for low and high initial density of cells is revealed in [21].

In this paper, a computational model based on the finite element method is proposed aiming at investigating the behavior of the two-dimensional chemotaxis–diffusion–convection system. We present numerical examples of different states of the system: (i) diffusion dominant, (ii) chemotaxis dominant, and (iii) convection dominant with the formation of descending plumes. Furthermore, we show that the chemotaxis sensitivity () and the taxis Rayleigh number () are two relevant parameters for the generation of instabilities. The effects of initial conditions and the initial cell density were also explored. Distinct initial settings lead to different solutions with plume patterns. Numerical tests show how the specified deterministic initial condition influences the overall behavior such as the number of plumes. In addition, the chemotaxis–diffusion–convection system is compared with the Rayleigh-Bénard and double-diffuse convection systems. The dimensionless parameters introduced in previous papers were renamed for better readability (table 2).

In section 2 of the present paper, the mathematical formulation of the chemotaxis–diffusion–convection system is presented. Section 3 gives the numerical method. Section 4 demonstrates the role of chemotaxis, in addition to a discussion on the simulation results and comparison with other models of buoyant convection. Finally, section 5 presents some concluding remarks.

2 Mathematical model

A mathematical model for the chemotaxis–diffusion–convection is proposed in [3] and reads as follows:


where denotes the velocity field of water (solvent), the pressure, and the water density and viscosity, the areal number density of bacteria (i.e., number of bacteria per unit area in a 2D space), and the volume and volumetric mass density of a bacterium, the concentration of oxygen, the buoyancy force exerted by a bacterium on the fluid in the vertical direction (unit vector ), the dimensional chemotaxis sensitivity, and the bacterium and oxygen diffusivities, the bacterial oxygen consumption rate, and the dimensionless cut-off function for oxygen concentration. The cut-off function is defined by the step function


where .

Bacteria are slightly denser than water and are diluted in the solvent, so that we consider and , respectively. The consumption of oxygen is proportional to the bacterial population density . Both and are advected by the fluid. When the oxygen concentration is lower than a threshold, the bacteria become quiescent [3], that is, they neither consume oxygen nor swim toward sites of higher oxygen concentrations. The dynamics of space filling, intercellular signaling, and quorum sensing are also neglected.

Figure 1: Boundary conditions for the system of equations (1-4). The air–water interface, where the oxygen concentration is equal to that of air, is not crossed by bacteria; the fluid vertical velocity component equals zero and the fluid is assumed to be free of tangential stress. The container walls are impermeable to bacteria and oxygen; a no-slip condition is imposed.

The system of equations (1-4) with the boundary conditions introduced in previous papers (e.g. [2, 18, 21]) is solved in a two-dimensional rectangular container (). The top boundary () represents the interface between liquid and air. On the free surface the concentration of oxygen is equal to the air concentration of oxygen () and the free tangential stress condition as well as the absence of bacterial flux are prescribed (figure 1). Therefore:


where is the unit outward normal vector. On the container walls (), a no-slip boundary condition is prescribed and the fluxes of bacteria and oxygen equal zero:


A no-slip condition at the air–water interface would enable the formation of hydrodynamic instabilities. The effect of a moving boundary due to the advection caused by an incompressible fluid flow is to be explored.

3 Computational model

3.1 Scaling and setting for numerical simulations

The characteristic length is defined by the container height and the characteristic bacterial density by the average of the initial bacterial population density defined as


This particular choice of characteristic bacterial density allows us to easily measure the total number of bacteria in each simulation for different initial distributions of bacteria.

Dimensionless variables are defined as


Five dimensionless parameters given below characterize the hydrodynamic and chemotaxis transport equations:


where is the taxis Prandtl number, the taxis Rayleigh number (buoyancy-driven flow), the taxis Lewis number, the dimensionless chemotaxis sensitivity, and the chemotaxis head. The taxis Prandtl, Rayleigh, and Lewis numbers are analogous to the respective heat and mass Prandtl, Rayleigh, and Lewis numbers in heat and mass transfer (table 2). The chemotaxis sensitivity () and head () characterize the chemotaxis system. Only and depend on the characteristic length and characteristic bacterial density .

After removing the prime in dimensionless quantities, the set of dimensionless equations becomes


The system (18) is solved in a rectangular domain with the initial conditions:


On the top of the domain , the dimensionless boundary conditions are prescribed as


while on the other boundaries, we impose the dimensionless boundary conditions as follows


In the following sections, we will refer to the hydrodynamic system for the first two equations in (18) and to the chemotaxis system for the last two equations in (18).

3.2 Numerical methods

The governing equations in (18) are solved using the finite element method. We adopt the biquadratic quadrilateral elements for the primitive variables and the bilinear quadrilateral elements for the primitive variables so as to satisfy the LBB (Ladyžhenskaya [22] - Babuška [23] - Brezzi [24]) stability condition. There are nine nodes in one biquadratic quadrilateral element and four nodes in one bilinear quadrilateral element. In each element, the function can be written as , where are the nodal unknowns.

To avoid the convective instability while solving the convection dominated flow equations, we adopt the idea of a streamline upwind/Petrov-Galerkin (SUPG) method [25]. We implemented an inconsistent Petrov-Galerkin weighted residual scheme developed in [26]. The resulting weak formulation reads as


In the context of inconsistent formulation, we have two kinds of test functions. For the non-convective terms we choose the test function to be the shape function . The test function is only applied to the convective term so as to add a numerical stabilizing term. The main purpose of employing a SUPG scheme is to add an amount of streamline artificial viscosity. The test function is therefore rewritten in two parts . is called the biased part. On each node , the biased part is where corresponds to the Cartesian coordinates and .

From the exact solution of the convection–diffusion equation in one dimension, following the derivation in [26], the constant is determined as


We define with being denoted as the grid sizes. Finally, the derived expression for is given below


For comparison purpose in section 4.4, the double diffusion system (26) and the Rayleigh–Bénard system (27) are solved using the software Freefem++ [27]. The code uses a finite element method based on the weak formulation of the problem. Taylor-Hood elements are chosen in FreeFem++. These elements are used together with a characteristic/Galerkin formulation to stabilize the convection terms.

3.3 Numerical validation of the coupled NS-KS equations

In this validation study, the following dimensionless differential equations accounting for the coupled Keller-Segel and incompressible viscous hydrodynamic equations are solved


The physical properties are set at the constant values of in . Equations in (18) are solved at in the continuously refined four meshes with . The predicted errors between the simulated and exact solutions, which are , , , and , are cast in their . The source terms , , and are derived from the previous exact solutions. From the predicted error norms, the spatial rates of convergence are plotted in Fig. 2.

Figure 2: The computed rates of convergence (roc) for the coupled set of Navier-Stokes and Keller-Segel equations. (a) roc= 3.95382 for , roc= 3.92594 for , roc= 1.96484 for ; (b) roc= 2.81384 for , roc= 1.80646 for .

The good agreement between the exact and simulated results and states of convergence demonstrate the applicability of the proposed SUPG scheme and the flow solver described in section 3.2 to investigate the chemotactic phenomena in hydrodynamic environment.

4 Numerical results and discussion

The linear stability analysis of the system (18) showed that the steady state becomes unstable for a range of physical parameters [2]. For sufficiently large characteristic bacterial density and Rayleigh number , hydrodynamic instabilities appear in the region near the surface at which the bacterial density is high. This instability may be related to the Rayleigh-Bénard instability occurring in thermal convection [3, 17]. This instability develops into a descending family of bacterium-rich plumes and leads possibly to the formation of convection cells (figure 3 (a)).

(a) CDC
(b) RBC
(c) DDC
Figure 3: Examples of convection cells for the chemotaxis–diffusion–convection (CDC), Rayleigh-Bénard convection (RBC), and double diffusive convection (DDC).

When is small, small perturbations of the velocity field are damped due to stabilizing effects of viscous friction and chemotaxis and therefore convective motion is negligible. It is shown that for low initial bacterium cell density the system (18) evolves to a steady and homogeneous state in the horizontal direction governed by the chemotaxis system [21, 3]. Like in the Keller-Segel or angiogenesis system, the random diffusion of cells in this case is balanced by the chemotaxis sensitivity of cells.

4.1 Descending plumes

In the present numerical simulations, the descending plumes can be described using three phases. In the first phase (figure 4 (a)-(c)), chemotaxis is a dominant mechanism. As bacteria consume oxygen, an oxygen concentration gradient is created that in turn provokes a bacterium motion toward the open surface where the oxygen is abundant. Bacteria chemotaxis causes the fluid to set in motion and counter rotating vortices to form (figure 3 (a)). Quickly, the bacterial density becomes quasi-homogeneous in the horizontal direction and is structured in layers in the vertical direction. A three-layer configuration is induced. A layer of higher concentration of bacteria forms below the air surface: the stack layer. When bacteria have migrated in the stack layer, a depletion layer is generated. At the container bottom, some bacteria are inactive, because the oxygen concentration decreases below a certain threshold. An inactive layer is established.

(a) t=0 (b) t=0.006 (c) t=0.012 (d) t=0.018 (e) t=0.024 (f) t=0.030 (g) t=0.036 (h) t=0.042
Figure 4: Evolution of the cell density at different times. Descending plumes of bacteria develop from an initial randomly distributed bacterium population.

The second phase (figure 4 (d)-(e)) exhibits high bacterial density on the surface. As bacteria swim toward air-supply region, the bacterial density in the stack layer increases. Consequently, advection from the counter rotating vortices becomes significant and brings in perturbations to the stack layer (figure 5). Therefore, advection causes instabilities to occur in the bacterial density distribution in the stack layer that sets the fluid in motion. At the same time, the oxygen concentration falls to its minimum at the bottom.

Figure 5: Evolution of the cell density number at the surface. In the initial stage of the chemotaxis–diffusion–convection, is homogeneous in the horizontal direction. As the cell density increases, hydrodynamic instability arises and the bacterial density becomes higher at some particular points.

In the third phase (figure 4 (f)-(h)), in fluid regions with a greater bacterial density in the stack layer, buoyancy force constrains bacteria to descend in the fluid. As a result, descending plumes of bacteria develop at these particular locations. This mechanism is analogous to the Rayleigh-Bénard instability in heat transfer problems where the fluid with a higher temperature and thus a lighter density than that above it rises to develop plumes of hot fluid.

4.2 Stabilizing effect of chemotaxis

When the average initial cell density is large, hydrodynamic instabilities appear in the system (18). This section is aimed at estimating and dimensionless governing parameters for hydrodynamic instabilities to appear as well as at analyzing the time scales for each of the three competitive physical mechanisms: (1) chemotaxis, (2) diffusion, and (3) convection of bacteria.

Convection is determined by the properties of the hydrodynamic system as well as the container height and via the taxis Rayleigh number . When the Rayleigh number is below the critical value for the hydrodynamic system of interest, motion of bacteria is primarily governed by diffusion and chemotaxis; when the Rayleigh number exceeds a critical value, bacterial taxis is primarily governed by convection [2].

Buoyancy force enables a fluid volume with a low bacterial density to ascend and a fluid volume with a high bacterial density to descend. Nevertheless, bacterium chemotaxis and friction dampen altogether the displacement.

The time scale for bacterium diffusion over the length scale is given by


The buoyancy force is balanced by friction in the fluid. Therefore, the time scale for the convective displacement of bacteria over the length scale can be defined by


The chemotaxis system is controlled by a competition between diffusion and chemotaxis of bacteria. Therefore, the chemotaxis time scale is expressed as


In a convection-dominant process, the convection time scale is smaller than the diffusion and chemotaxis time scales (table 3). Hence:


In particular, . The value of the critical Rayleigh number associated with the convection given by [2] in their linear stability analysis is described by the solution of an ordinary differential system. The second condition on the time scales leads to the following inequality


Similarly to the critical Rayleigh number (), a critical number can be introduced.

On the other hand, the chemotactic motion is predominant when the chemotaxis time scale is smaller than the diffusive time scale and convection is negligible:

Dominant Dominant Dominant
convection diffusion aerotaxis
and ;
Table 3: Phenomenological analysis based on time scales of the three competitive mechanisms: chemotaxis, diffusion, and convection of bacteria. Both oxygen and bacterial diffusion favor fluid homogenization.

The effects of and the product were tested subject to the random initial condition. The results are shown in figure 6 for and . Our numerical results agree with the linear stability analysis carried out in [2]. at first falls as increases to reach its minimum value and then rises again. A sufficient decrease or increase of may promote stabilization as both taxis and diffusion operate as the fluid-homogeneization factors, either directly () or indirectly (). When is small, stabilizing effect is due to bacterial diffusion and the solution is of the diffusive type. When is large, stabilizing effect results from the bacterial taxis and the solution is of the chemotactic type.

Figure 6: The stable and unstable regions of the system (18) are plotted in terms of the critical taxis Rayleigh number and the product of the dimensionless chemotaxis sensitivity and chemotaxis head . The points correspond to the predicted values of . The red line corresponds to a fit of data in two parts: on the left side it is fitted by a power function and on the right side by a linear function. In the unstable region, the predicted solution is of the convective type (b). In the stable region the predicted solution can be of the diffusive type (a) or the chemotactic type (c).

The taxis Rayleigh number () that characterizes the competition between diffusion and convection plays the same role as the Rayleigh number plays in classical convection. When the Rayleigh number increases, the gravitational force becomes predominant. When rises, competition between chemotaxis and convection of bacteria becomes stronger. The condition (23) suggests that increases linearly with respect to , as illustrated in figure 6. We can affirm that chemotaxis has a stabilizing effect on the differential system of current interest.

4.3 Distribution and number of plumes and initial conditions

4.3.1 Position and spacing of plumes

The exact localization of plume generation was investigated. Under the random initial condition introduced in [21], the location is hard to be predicted. We consider deterministic initial conditions to investigate the formation and location of descending plumes. Many initial conditions among the set of tested distribution of bacterial settings are able to trigger convective patterns.

All numerical results are computed at , , , , and to ensure the formation of bacterium rich plumes. We first consider a profile given by the form of wave function to determine the two initial layers of bacteria. Simulation results and initial conditions are shown in figure 7. The upper layer has a higher bacterial density than the lower layer. The plumes form at the initial location of the crest of the wave function where there is a larger number of bacteria. By reversing the initial conditions with a denser inferior layer (figure 8), the plumes also form at the same location. We can observe that distinct initial conditions can lead to the formation of a very similar pattern (figures 7-8).

Figure 7: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the upper layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).
Figure 8: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the upper layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).

We vary the profile of the curve between the upper and lower layers in the initial condition with different wave functions , , and . The wavenumber is defined as the spatial frequency of waves per unit distance. The wavelength is defined as the domain horizontal length divided by the wavenumber. Increasing the wavenumber of the initial wave function causes the formation of additional plumes to occur (figures 9 and 10). In all simualtions, the locations of the plumes correspond to sites of initially greater local bacterial density.

Figure 9: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the upper layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).
Figure 10: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the lower layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).

When the wavenumber in the initial profile is large enough, the number of plumes is not equal to the wavenumber fixed by the initial condition. Three plumes, at most, form (figures 11 to 14). Nonetheless, more than three plumes can form, but they merge later (figure 13). Plume merging was previously observed in [21]. The plume merging mechanism is analogous to that of the Rayleigh-Bénard flow. Spacing between plumes seems to be intrinsic to the system and the position can only be predicted in very specific cases such as the examples given above.

Figure 11: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the upper layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).
Figure 12: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the lower layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).
Figure 13: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the upper layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).
Figure 14: Numerical results for with respect to the deterministic initial conditions. At the initial time, bacterial density is higher in the lower layer. (a) initial condition is ; (b) simulation result at for the initial condition given in (a) ; (c) initial condition is ; (d) simulation result at for the initial condition given in (c).

The randomly perturbed bacterial density is defined as follows [21]:


with being a random number with a uniform probability distribution over . Simulation results with the randomly perturbed initial condition are given in figure 15. We note that the plume-to-plume spacing is not fixed but the number of plumes in the solution remains the same as that with the solution obtained subject to the deterministic initial condition illustrated in figures 11 to 14.

Figure 15: Simulations subject to the initial condition given in (25) (left) and the corresponding numerical results for at (right).

In figures 715, we can see plumes forming at the border of the domain. These border plumes seem to be caused by the no-slip boundary condition on the velocity depending on the arrangement of the convection cells. Bacteria first agglomerate at the surface near the wall. Because the fluid velocity is equal to zero at the wall boundary, a large amount of bacteria remains close to the wall. At a location away from the wall, the velocity is non-zero, thus bacteria descend into the fluid and form a plume at the border.

The location of the plumes may only be predicted in very simple tests for which the wavelengths of initial conditions are lower than the wavelength of the process. Despite the difficulty to predict the site of plume formation, the system presents a dominant wavelength of the fingering instability. The wavelength and wavenumber are redefined such as those given in [28]. The wavelength is the length of the domain divided by the number of plumes and the wavenumber divided by the wavelength. The wavelength and wavenumber are similar for all simulation tests whatever the length of the domain is (tables 4 to 7).

number of plumes wavelength wavenumber
2 2 3.14
2 2 3.14
2 2 3.14
2 2 3.14
Table 4: The predicted number of plumes, wavenumber and wavelength for . Each row corresponds to a different simulation result subject to the randomly perturbed initial condition given in (25).
number of plumes wavelength wavenumber
3 2 3.14
3 2 3.14
4 1.5 4.19
3 2 3.14
Table 5: The predicted number of plumes, wavenumber and wavelength for . Each simulation result is subject to the random initial condition given in (25).
number of plumes wavelength wavenumber
5 1.6 3.93
5 1.6 3.93
5 1.6 3.93
5 1.6 3.93
Table 6: The predicted number of plumes, wavenumber and wavelength for . Each simulation result is subject to the random initial condition given in (25).
number of plumes wavelength wavenumber
5 2 3.14
6 1.67 3.77
5 2 3.14
5 2 3.14
Table 7: The predicted number of plumes, wavenumber and wavelength for . Each simulation result is subject to the random initial condition given in (25).

4.3.2 Growth of plumes

We now arbitrarily define the plumes by the isoline . This choice allows us to track the stack layer and the plumes being formed. The layer where represents the layer with a high bacterial density from which plumes form and descend in the fluid. On the other hand, the layer with corresponds to the depletion layer, from which bacteria have migrated to the stack layer.

We define the growth rate of the plume amplitudes by . The plume amplitude is computed by measuring the distance from the surface to the tip of the plume at time (table 8) and is the time increment. Data obtained can be interpolated by the exponential function as the amplitudes of the descending plumes grow exponentially (figure 16). The exponential growth is well understood in the conventional Rayleigh-Taylor convection. Subsequent to the exponential growth phase, the growth rate of amplitude decreases, as plumes get closer to the bottom of the container.

Figure 16: Plot of the growth rate of plume amplitude (dots) corresponding to the data in table 8. The growth rate is interpolated by the exponential function (line).
(a) (b)
0.18 0.16 0.18 0.16
0.19 0.18 0.46 0.20 0.18 1.07
0.20 0.19 0.76 0.21 0.19 1.07
0.21 0.20 0.92 0.23 0.20 1.99
0.23 0.21 1.53 0.28 0.21 3.53
0.26 0.22 2.91 0.37 0.22 8.13
0.34 0.24 6.44 0.57 0.24 16.34
0.51 0.25 13.81 0.78 0.25 17.12
0.72 0.26 17.81 0.87 0.26 8.13
0.83 0.27 9.67 0.93 0.27 4.91
Table 8: The predicted growth rate of the plume amplitudes for two simulations subject to the random initial condition (25).

4.4 Comparison with other buoyancy-driven convections

The chemotaxis–diffusion–convection system has many features similar to other well known buoyancy-driven flows, such as the double diffusive and Rayleigh Bénard convection.

Double diffusive convection occurs in a fluid containing at least two components with different diffusivities. A destabilizing component diffuses faster than the stabilizing one [29]. The distinct diffusivities yield a density difference capable of driving the motion of fluid [30].

Comparison with a classical example of double diffusive phenomenon in oceanography can be carried out by considering two superposed fluid layers with a specific combination of temperature () and the solute concentration (), namely the salinity. The system of dimensionless equations of the double diffuse problem is the following:


where and are the thermal and mass Rayleigh number, respectively, the Prandtl number, and the Lewis number (table 2).

The stability of the system that exhibits a diffusive and a finger regime depends on both Rayleigh number types. In the finger regime, a small perturbation at the interface between the layers develops into a pattern of descending fingers, for instance salt fingers.

Another buoyancy-driven convection is the Rayleigh-Bénard convection that arises by fluid in a reservoir heated from below. Convection results from thermal gradient. The set of equations is as follows:


and are the Rayleigh and Prandtl number, respectively. The balance between the gravitational and viscous forces is expressed by the Rayleigh number . When is larger than a critical value that can be obtained analytically, convective patterns appear [31].

The double diffusive convection equations in (26) and the Rayleigh-Bénard convection equations in(27) are solved using our numerical method (figures 17 and 18). Plumes formed in the Rayleigh-Bénard convection process are similar in shape to plumes of the chemotaxis–diffusion–convection. However, double diffusive and Rayleigh-Bénard convections exhibit both ascending and descending plumes, whereas chemotaxis–diffusion–convection is only characterized by descending plumes.

Figure 17: Plot of the solution in the double diffusive system (26). Red area designates the region rich in , while yellow area is poor in . Descending plumes rich in and ascending plumes poor in form at the interface of two layers of fluid. The green color designates a region where the concentration of is slightly higher than that in the yellow area due to diffusion. (a) Initial condition: , if , and , if . (b) at time =0.1 for , , and .
Figure 18: Plot of the solution in the Rayleigh-Bénard system (27). Descending plumes of lower temperature and ascending plumes of higher temperature form in the fluid. , .

In each problem, the effective Rayleigh numbers depend on the temperature difference between the opposite fluid domain surfaces, the gradients of temperature and salinity between the fluid layers, and the difference of density between the bacteria and solvent, for the Rayleigh-Bénard, double diffusive, and chemotaxis–diffusion–convection systems, respectively. All of the three diffusion–convection processes have similar and distinct features (table 9). The involved dimensionless parameters are listed in table 2.

I.C. for fingers Layers Any Layers
B.C. for mass Neumann Neumann Dirichlet
Physical regime Diffusive & Diffusive, Diffusive &
convective convective & convective
Table 9: Recapitulative of the physical mechanisms involved in the double diffusive convection (DDC), chemotaxis–diffusion–convection (CDC) and the Rayleigh-Bénard convection (RBC).

Figure 3 shows that the arrangement of the convection cell structure can be the same for the three systems mentioned above. A particular position of the domain will be in a clockwise cell in some simulations and counter-clockwise cell in others. This behavior is also observed in Rayleigh-Bénard convection simulations.

In the chemotaxis–diffusion–convection system, chemotaxis plays an essential role in the early stage, as it organizes the fluid domain. In a rectangular domain, starting from a given initial condition, aerotaxis of bacteria builds quasi-homogeneous layers in the horizontal direction, similarly to the double diffusive [30]. The multi-layered fluid creates a density gradient between the top of the stack layer and the bottom of the depletion layer that is similar to the temperature gradient set by the Dirichlet boundary condition imposed in the Rayleigh-Bénard case. The chemotaxis–diffusion–convection system evolves itself to a proper condition that leads to instabilities.

Chemotaxis, that is not present in the other systems, brings flexibility in the choice of initial as well as boundary conditions. An analogy between the differential systems should stimulate further mathematical analysis of the system and a better understanding of the role of chemotaxis in convection problems.

5 Concluding remarks

In this paper, we have studied the chemotaxis–diffusion–convection system with the focus on the differential system rather than on the experimental settings. Our studies are based on the numerical implementation of the upwind finite element method with inconsistent Petrov-Galerkin weighted residual scheme to solve the coupled convective chemotaxis-fluid equations. Several simulation examples have been presented. They exhibit physically different spatially organized convection patterns.

The chemotaxis–diffusion–convection system can be described by three regimes. In the convective regime, the formation of plume patterns proceeds a temporal sequence of development stages. First, a spatial layered structure is created. Bacteria agglomerate in the upper stack layer, thus forming a depletion layer, where the bacterial density is very low. Subsequently, bacterial convection strengthens with time and instabilities in the stack layer appear. Finally, plumes of bacterial falling in the fluid are generated. The growth rate of the plumes is of the exponential type. The wavelength of the growing instabilities is defined from the parameters of the problem. Initial conditions appear to have a small influence on the plume number. The location of plumes can only be predicted when a very simple initial condition is set up.

In the diffusive and chemotactic regimes, the chemotaxis system has a stabilizing effect on the fluid. When the taxis Rayleigh number increases, the gravitational force becomes dominant and instabilities occur and grow. From the phenomenological analysis and the numerical results, the critical Rayleigh number increases proportionally with the product of the chemotaxis sensitivity and head, when the chemotaxis sensitivity is large.

In all buoyancy-driven flows (chemotaxis–diffusion, double diffusive, and Rayleigh-Bénard convection), the dimensionless differential systems and the convection patterns are similar. Analogy between these types of convection should launch further analytical studies of the chemotaxis–diffusion–convection problem.


Y.D. received a PhD grand from Pierre and Marie Curie University - Sorbonne University and was also supported by the CASTS at National Taiwan University.


  • [1] M. P. Brenner, L. S. Levitov, E. O. Budrene, Physical mechanisms for chemotactic pattern formation by bacteria, Biophysical Journal 74 (4) (1998) 1677–1693.
  • [2] A. J. Hillesdon, T. J. Pedley, Bioconvection in suspensions of oxytactic bacteria: linear theory, Journal of Fluid Mechanics 324 (10) (1996) 223–259.
  • [3] A. J. Hillesdon, T. J. Pedley, J. O. Kessler, The development of concentration gradients in a suspension of chemotactic bacteria, Bulletin of Mathematical Biology 57 (2) (1995) 299–344.
  • [4] N. Hill, T. Pedley, Bioconvection, Fluid Dynamics Research 37 (1–2) (2005) 1–20.
  • [5] E. F. Keller, L. A. Segel, Model for chemotaxis, Journal of Theoretical Biology 30 (2) (1971) 225–234.
  • [6] C. S. Patlak, Random walk with persistence and external bias, The Bulletin of Mathematical Biophysics 15 (3) (1953) 311–338.
  • [7] A. Blanchet, J. Dolbeault, B. Perthame, Two-dimensional Keller-Segel model: optimal critical mass and qualitative properties of the solutions, Electronic Journal of Differential Equations 44 (2006) 1–33.
  • [8] M. A. Herrero, J. J. L. Velázquez, Singularity patterns in a chemotaxis model, Mathematische Annalen 306 (1) (1996) 583–623.
  • [9] W. Jager, S. Luckhaus, On explosions of solutions to a system of partial differential equations modelling chemotaxis, Transactions of the American Mathematical Society 329 (2) (1992) 819–824.
  • [10] T. Nagai, T. Senba, Global existence and blow-up of radial solutions to a parabolic-elliptic system of chemotaxis, Advances in Mathematical Sciences and Applications 8 (1998) 145–156.
  • [11] J. J. L. Velázquez, Point dynamics in a singular limit of the Keller-Segel model 2: Formation of the concentration regions, SIAM Journal on Applied Mathematics 64 (4) (2004) 1224–1248.
  • [12] B. Perthame, Cell motion and chemotaxis, in: Transport Equations in Biology, Birkhauser, Basel, 2007, pp. 111–149.
  • [13] M. Chaplain, S. McDougall, A. Anderson, Mathematical modeling of tumor-induced angiogenesis, Annual Review of Biomedical Engineering 8 (1) (2006) 233–257.
  • [14] N. V. Mantzaris, S. Webb, H. G. Othmer, Mathematical modeling of tumor-induced angiogenesis, Journal of Mathematical Biology 49 (2) (2004) 111–187.
  • [15] L. Corrias, B. Perthame, H. Zaag, A chemotaxis model motivated by angiogenesis, Comptes Rendus Mathematique 336 (2) (2003) 141–146.
  • [16] L. Corrias, B. Perthame, H. Zaag, Global solutions of some chemotaxis and angiogenesis systems in high space dimensions, Milan Journal of Mathematics 72 (1) (2004) 1–28.
  • [17] A. Metcalfe, T. Pedley, Bacterial bioconvection: weakly nonlinear theory for pattern selection, Journal of Fluid Mechanics 370 (1998) 249–270.
  • [18] I. Tuval, L. Cisneros, C. Dombrowski, C. W. Wolgemuth, J. O. Kessler, R. E. Goldstein, Bacterial swimming and oxygen transport near contact lines, Proceedings of the National Academy of Sciences of the United States of America 102 (7) (2005) 2277–2282.
  • [19] R. Duan, A. Lorz, P. Markowich, Global solutions to the coupled chemotaxis-fluid equations, Communications in Partial Differential Equations 35 (9) (2010) 1635–1673.
  • [20] J.-G. Liu, A. Lorz, A coupled chemotaxis-fluid model: Global existence, Annales de l’Institut Henri Poincare (C) Non Linear Analysis 28 (5) (2011) 643–652.
  • [21] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, P. A. Markowich, Sinking, merging and stationary plumes in a coupled chemotaxis-fluid model: A high-resolution numerical approach, Journal of Fluid Mechanics 694 (2012) 155–190.
  • [22] O. A. Ladyženskaja, The mathematical theory of viscous incompressible flow., Gordon and Breach Science, New York, 1969.
  • [23] P. I. Babuška, Error-bounds for finite element method, Numerische Mathematik 16 (4) (1971) 322–333.
  • [24] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers, ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique 8 (R2) (1974) 129–151.
  • [25] A. N. Brooks, T. J. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 32 (1–3) (1982) 199–259.
  • [26] M.-T. Wang, Development of finite element method for incompressible Navier-Stokes equations, Ph.D. thesis, National Taiwan University (Jul. 1996).
  • [27] F. Hecht, New development in FreeFem++, Journal of Numerical Mathematics 20 (3-4) (2013) 251.
  • [28] T. W. Pan, D. D. Joseph, R. Glowinski, Modelling Rayleigh–Taylor instability of a sedimenting suspension of several thousand circular particles in a direct numerical simulation, Journal of Fluid Mechanics 434 (2001) 23–37.
  • [29] L. Lemaigre, M. A. Budroni, L. A. Riolfo, P. Grosfils, A. D. Wit, Asymmetric Rayleigh-Taylor and double-diffusive fingers in reactive systems, Physics of Fluids 25 (1) (2013) 014103.
  • [30] J. S. Turner, Double-diffusive phenomena, Annual Review of Fluid Mechanics 6 (1) (1974) 37–54.
  • [31] E. Bodenschatz, W. Pesch, G. Ahlers, Recent developments in Rayleigh-Bénard convection, Annual Review of Fluid Mechanics 32 (1) (2000) 709–778.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description