Team:KU Leuven/Project/Ecological/Modelling

Wind speed

Because of friction and obstacles on the earth’s surface, wind speed varies with altitude. Generally, the velocity decreases with increasing altitude. A logarithmic wind profile is appropriate for the part above the crops (Goudriaan, 1977, p. 96). The formula for this profile is

with u representing the velocity. Here d accounts for an upward shift above a vegetative cover. The relation d=0.63 x z_c is suggested, where z_c is the height of the crops. The length z₀ is called the roughness length and is often supposed to be about one tenth of z_c.

For the part inside the canopy, the profile is exponential.

As can be noticed, the wind speed decreases exponentially with extinction factor a when approaching ground level and thus going deeper into the canopy. u_c is the speed at height z_c and can easily be calculated by the formula of the logarithmic wind profile.

Evidently, the wind direction changes in time. Most regions however tend to have a dominant wind (“Prevailing winds,” Wikipedia, 2013), at least during a certain time period. This is the reason our model only incorporates convection in one direction, which greatly simplifies calculation. Furthermore, average wind speeds are measured and published monthly by the Royal Meteorological Institute (KMI) in Belgium (“Maandelijkse normalen - KMI,” 2013). These are measured at 10 m above the ground, making it possible to calculate u*/k used in the log law. Now, all parameters for convection are known when the crop height is given. The wind profile was entered in the software using a piecewise function and self-defined parameters, making it easy to change wind speed and crop height. An example of the profile is plotted in Figure 1.

Boundary conditions

In order to solve our model, boundary conditions need to be specified. We'll represent the environment as a square box (5X5X5m) with 6 faces. Within this environment box, we assume a bacterium box (5X5X1cm), suspended in "mid air", with its own 6 (smaller) faces.

Inflow and outflow

For the environment box, the faces perpendicular to the wind direction were specified as an inflow (with a concentration of 0 mol/m³), respectively an outflow.

No flux

Within the environment box, the face representing the ground will have no flux through it, because we assume that the diffusion coefficient of our pheromones is much larger in soil than in air. Within the bacterium box, all faces but one (pheromone outlet) are given zero flux as well.

Open boundaries

The remaining 3 faces of the environment box remain "open" since the pheromone flux through these faces is unknown at this point.

Specified flux

The upper face of the bacterium container is the only place where pheromones are released in the air. Under the assumption that the pheromone production and vaporisation is at steady state, we can set the flux at this surface equal to the production rate of the entire colony.

Once we know the production per cell (by the procedure in MeS Modelling) it is straightforward to calculate the output of the whole colony by using an average cell density in the appropriate growth phase. This output represented as an amount per time can then be converted to a flux through the contact surface with the air, which can be entered in the simulation program.

Results

Although we lacked real values for the production rates of the pheromones, we ran the model with a "ballparked" flux of 1 mmol/(m²s) to get an idea of the shape of the solution. We used a simple raised cosine function to capture the oscillating behaviour of the oscillator with a period of 1 hour. As can be seen in the animated figure below, convection by wind is the dominant effect on dispersion of the pheromones. Therefore the solution follows the oscillations quite well and lag or local accumulation of pheromones is not a problem. Note that the solution has negative values at some points, which is physically impossible of course. This is solely due to numerical errors and can be avoided by the use of a finer mesh (at a much larger computational cost).

Software

At first, we used Mathworks Matlab with the Partial Differential Equations Toolbox, but this software was limited to 2D geometries and more stringent boundary conditions. Later on, we noticed that the Department of Chemical Engineering of our university provided COMSOL Multiphysics. This program was very suitable for our purposes, as it provides a model for “Transport of Diluted Species”. Our complete model can be downloaded here.