May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Particle Tracking Model for Simulation of Advection and Diffusion in Prairie Dog Burrows

From MathBio

Jump to: navigation, search

Executive Summary

This project uses a stochastic Lagrangian simulation to approximate the advection and diffusion behavior in animal burrows. This process is governed by a second-order partial differential equation, which is not always easy to solve. Instead, individual molecules are programmed to move randomly, and the overall trends are used the draw conclusions. These random paths are simplified using well-understood statistical approximations for Brownian motion. The molecular model is supplemented with a fluid mechanics model to account for advection.

This model allows for the study of molecular diffusion and advection, and how these affect animal respiration. For an animal to survive, oxygen must continuously flow in and carbon dioxide must continuously flow out. With diffusion alone, the author concludes that large animals would be unable to survive in underground burrows.


The exchange of gases within an animal burrow has been an extensive area of study due to the barrier that the burrowing medium (both aquatic and terrestrial) often imposes to the diffusion of gases. While smaller animals run into limited problems with meeting their necessary aerobic metabolic needs through gaseous diffusion, larger animals are often unable to survive on diffusion alone. For example, for a 1 meter deep burrow the average diffusion time of molecular oxygen approaches 27 hours, leading to a flux of 3\times 10^{{22}} oxygen molecules per hour, two orders of magnitude less than is needed for respiration. In addition, the buildup of hazardous substances, byproducts of respiration, have just as much difficulty escaping the burrow after being released. Studies of mammalian burrows have observed carbon dioxide concentrations exceeding 15 percent, extremely high when compared with normal atmospheric levels of only .0387 percent. Therefore, alternative mechanisms of gas exchange are needed in order for these large animals to survive within their burrows. Advective flushing is one such mechanism used to mediate gaseous exchange, a process much more rapid than gaseous diffusion. Prairie dog burrows are often cited as the gold standard as they take full advantage of advection due to their means of construction, which will be the focus of this review.

Rather than utilizing analytical solutions to the advection-diffusion equation,

{\frac  {\partial C}{\partial t}}+u{\frac  {\partial C}{\partial x}}-{\frac  {\partial ^{2}C}{\partial x^{2}}}=0

to compute concentrations within the burrow, a stochastic particle dispersion model will be used to simulate the dispersion of carbon dioxide within the burrow. By releasing a large number of particles, and applying the basic physics of dispersion and advection to each particle, mean concentrations can be computed that match closely with analytical solutions. The basics of this model will be discussed below, followed by the results of these simulations.

Simulating the Advection-Diffusion Equation

In 1826, the botanist Robert Brown noticed random fluctuations in the pollen grains that he was studying under the microscope. This erratic behavior at the microscopic scales would later become known as Brownian motion. At the time, Brown believed that he had found a universal force that governs movements of living organisms. It wasn't until the discovery of the atom that this behavior could be accounted for by the random bombardment of molecules against the pollen grains. In 1906, Albert Einstein submitted a paper that quantified this motion. Einstein was interested in how far a particle would travel in a given time interval. Due to the large number of collisions, classical physics was unable to describe this behavior. Using the diffusion equation in combination with a moments analysis, however, led Einstein to discover that the mean square displacement of a particle could be expressed as,

{\bar  {x^{2}}}=2Dt

where {\bar  {x^{2}}} is the mean squared displacement (second moment), D is the diffusivity of the substance being considered and t is the amount of time elapsed. Knowing this relation, scientists could now estimate the growth of a cloud of particles due to molecular diffusion.

Numerical simulation of the advection-diffusion equation was first proposed by the fluid dynamisist G.I. Taylor. In his 1921 paper, "Diffusion by Continuous Movements", Taylor acknowledged that diffusion at a large scale (such as the atmosphere) could be modeled in the same way as the molecular scale. Taylor went so far as to suggest one of the first Monte Carlo techniques when he suggested,

"Suppose that a point starts moving with uniform velocity v along a line, and that after a time T it suddenly makes a fresh start and either continues moving forward with velocity v or reverses its direction and moves back over the same path with the same velocity v. Suppose that this process is repeated n times and that we consider the mean values of the quantity concerned for a very large number of such paths."

An example of this reasoning is shown in the figure below. As more particles are released, the plume fills in and convergence can begin. One this plume fills in with hundreds of thousands of particles, it would be clear that the length scale of this plume grows as {\sqrt  {2Dt}}.


As computer power advanced during World War II, the ability to simulate the Advection-Diffusion equation using a large number of independent particles increased dramatically. The motion of an individual particle is given by

\Delta x=u\Delta t+{\sqrt  {2Dt}}~N(0,1)

where N(0,1) is a random number drawn from a normal distribution with a mean of 0 and a variance of 1, and u is the advective velocity, which acts as a memory term. This movement was calculated at discrete time steps (\Delta t) for as many as 1 million particles. To calculate concentrations, a grid was overlaid onto the particle field, and the number of particles within each grid box was then counted. Each particle would be given a predetermined mass, which was determined by the amount of material being diffused. By summing the particles in each grid box, multiplying by the mass of each particle, and then dividing by the volume of the grid, concentrations at a given location could be calculated. Particle simulations used within this study are shown below. The top plots represent pure diffusion, while the bottom plots include advection and diffusion within the .1 m wide burrow. Time advances to the right, t(3) > t(2) > t(1).


Model Formulation and Results

For the prairie dog burrow, the following parameters were chosen for the particle diffusion model:

  1. The length of the burrow L = 10 m
  2. The width of the burrow W = .1 m
  3. CO_2 respiration rate for a .5 kg mammal = 100 ml/min.

This extended study focused on three different scenarios: 1) diffusion through a burrow with impermeable walls 2) diffusion through a burrow with permeable walls and 3) Advection and diffusion through a burrow with impermeable walls. For each scenario, the initial particle distribution was a line source extending across the width of the burrow. Since the prairie dog fills much of it's burrow, this well mixed criterion is justified.

Diffusion with impermeable walls

To simulate the continuous release of carbon dioxide, a line source beginning at x=0 was released after every time step, as seen in the particle figure above. The particles within this line source were then transported according to the numerical equation discussed earlier. When a particle encountered a wall, it was reflected from the boundary. This simulation was performed over the course of 1 hour, and is shown below. From this figure, it is clear that the build up of carbon dioxide would lead to rapid suffocation, even after only 1 hour of respiration in the burrow. The concentration near the prairie dog approaches 3\times 10^{3} mg/m^3, which is equivalent to 2\times 10^{5} ppm, or 20 percent of the mixture! This further reinforces Wither's conclusion that large animals would be completely unable to survive on diffusion alone.


Diffusion through porous media

To include porous media, the diffusion coefficient had to be altered when a particle transitioned between the air of the burrow and the soil. From the review of the Withers' article, we saw that the altered diffusion coefficient was a function of the soil porosity (f), which ranges between 0 and 1 and is given by,

f={\frac  {V_{a}+V_{w}}{V_{a}+V_{w}+V_{s}}}

where V_{i} is the volume of air, water, and soil respectively. See MBW: Gas Exchange in Animal Burrows for more detail on this calculation. The altered diffusion coefficient is then given by,


In addition to altering the diffusion coefficient, I also had to create an algorithm that would determine whether a particle reflected off of the burrow walls or was absorbed into the soil. Since the soil porosity parameter f ranges between 0 (impermeable) and 1 (completely permeable), this can be used as a probability that a particle is reflected or absorbed. Whenever a particle contacted a wall, a random number was drawn between 0 and 1. If that random number was less then f, the particle was absorbed into the soil.

For this part of the study, I looked at the effect that soil porosity had on the buildup of carbon dioxide in the burrow, which is shown in the figure below. From this figure, it is clear that soil porosity can greatly enhance the ventilation of the burrow. As the soil becomes more porous (f increases), carbon dioxide has more routes to diffuse out of the burrow, leading to lower concentrations within the burrow. Even for completely permeable soils, however, the concentrations are still much too large for prairie dog survival, as the carbon dioxide percentage approaches 5 percent of the mixture.


Advection and Diffusion

The implementation of the advection-diffusion model was done in the same way as the pure diffusion model, but with the inclusion of the advective u\Delta t memory term discussed earlier. Three advective velocities were chosen to simulate the ventilation of carbon dioxide (.001 m/s, .01 m/s, and .1 m/s). Studies have confirmed that velocities within the burrow of prairie dogs reach a maximum of 1 m/s, making these velocity scales relevant estimations. In order to obey the no-slip condition of fluid mechanics, a parabolic velocity profile was chosen that the velocity goes to 0 at the boundaries and is maximum (u_{0}) along the centerline of the burrow,

u(y)=4u_{0}{\frac  {y}{W}}{\bigg (}1-{\frac  {y}{W}}{\bigg )}

where W is the width of the burrow and y is the position of the particle. The figures of these three advective scenarios are shown below. For the u_{0}=.001 m/s case, concentrations are reduced to approximately 50,000 mg/m^3, a mixing ratio of 2.5 percent. Studies on humans have shown that any carbon dioxide level above 1 percent will result in unconsciousness. Therefore, these prairie dogs will be unable to survive with an advective velocity below .001 m/s. For the other two scenarios, concentrations of carbon dioxide are reduced to survivable levels of 5,000 mg/m^3 (.2 percent) for a velocity of .01 m/s and 500 mg/m^3 (.03 percent) for a velocity of .1 m/s. These carbon dioxide levels are roughly what would be measured within an office building and the atmosphere, respectively. The significance of advective flushing can be gathered by observing all three plots and noticing how each case reaches a steady state with an advancing front that progress along the burrow with increasing time. Because of advection, concentrations do not build in the burrow as they do for pure diffusion.

Adv1.png Adv2.png Adv3.png

Temperature Effects

Ganot et al. expanded the work on convection-advection in animal burrows to include temperature effects. During the morning and evening, the temperature gradient across the burrow is sufficient to induce Thermal Convective Venting (TCV.) Previous studies have measured this air flow to be in the range of 0.2 m/s.

To test these effects, the authors used a combined experimental-theoretical approach, choosing Meriones crassus as the case study mammal. The temperature profile found experimentally was inconsistent with a system dominated by conduction alone, so the authors conclude that convection must take place as well. To mathematically model the system, the authors employed a numerical solution to the heat equation in axisymmetric cylindrical coordinates.


Laboratory measurements of TCV show that the convection effects are more than sufficient to account for the carbon dioxide production of a small mammal, even without considering temperature generation from the animal’s metabolism. Mammals larger than 10 kilograms might have sufficient fluxes during upset conditions.


From this study, we have concluded that large animals are unable to survive on diffusion alone, as carbon dioxide concentrations will rise to levels that result in death, with mixing ratios ranging between 5 and 20 percent. By adding advection through the burrow, however, concentrations are reduced significantly, allowing for ventilation that is sufficient enough to keep carbon dioxide levels below the danger threshold, as long as the advective velocity is greater than .01~m/s. The necessity of active ventilation can be observed in nature, as prairie dogs construct their mounds in such a way as to promote advective flushing.


Ganot, Y., M. I. Dragila, and N. Weisbrod. "Impact of thermal convection on air circulation in a mammalian burrow under arid conditions." Journal of Arid Environments 84 (2012): 51-62.

Withers, P.C. "Models of Diffusion-Mediated Gas Exchange in Animal Burrows," The American Naturalist. vol 112, pp. 1101-1112, 1978

Matlab Code For Particle Diffusion Model

clear all
close all

T = 3600; %1 hour in burrow
dt = .1;

Np = 500; %particles per timestep
W = .1;  %Width of Burrow
L = 10; %Length of Burrow

%Create a grid to count particles

dx = .01;
dy = .01;

[x,y] = meshgrid(-L/2:dx:L/2,0:dy:W);

u_0 = -.1; %Advective Velocity

D = 16E-6; %Diffusivity of C02 in air at 20 deg C

rho_C02 = 1.98E6; %Density of C02 mg/m^3

VC02 = (100/(1000*1000*60))*rho_C02;   % typical respiration for this size mammal is 100 mL/min, which then is converted to mg/s

masspart = (VC02*dt)/(Np);

xpartini(1:Np) = 0;
ypartini = linspace(0.0,W,Np);

xpart = xpartini;
ypart = ypartini;

h = waitbar(0,'Please wait...');

i = 0;

%particle counting algorithm

bin(1:size(x,1),1:size(y,2)) = 0;

for t = dt:dt:T
    i = i + 1;
u = 4.*u_0.*(ypart./W).*(1-ypart./W); %Advective Velocity;

xpart = xpart + u.*dt + sqrt(2.*dt.*D).*randn(1,i*Np);
ypart = ypart + sqrt(2.*dt.*D).*randn(1,i*Np);

%Particle Reflection off of walls


xpart = [xpartini,xpart];
ypart = [ypartini,ypart];

%Uncomment the lines below to see simulation
% plot(xpart,ypart,'.','MarkerSize',3)
% xlim([-1 1])
% pause(0.01)

for count1 = 1:size(xpart,2)
   if xpart(1,count1) > x(1,size(x,2)) || xpart(1,count1) < x(1,1) 
    m = (size(x,2)-1)/(x(1,size(x,2))-x(1,1));
    k = m.*(xpart(1,count1)-x(1,1))+ 1;
    k = round(k);

    m = (size(y,1)-1)/(y(size(y,1),1)-y(1,1));
    j = m.*(ypart(1,count1)-y(1,1))+ 1;
    j = round(j);
    bin(j,k) = bin(j,k) + masspart;

concentration(:,:,i) = bin./(dx*dy);


x1 = -L/2:dx:L/2;

fh = figure(2);


hXLabel = ylabel('Concentration (mg/m^3) ');
hYLabel = xlabel('Burrow Length (m)');
titlabel = title('');

set([hXLabel, hYLabel,titlabel], ...
    'FontName'   , 'AvantGarde');
set([hXLabel, hYLabel,titlabel]  , ...
    'FontSize'   , 12         );

set(gca,'FontName' ,'Helvetica','FontSize',12, ...
  'Box'         , 'off'     , ...
  'TickDir'     , 'out'     , ...
  'TickLength'  , [.02 .02] , ...
  'XMinorTick'  , 'off'      , ...
  'YMinorTick'  , 'off'      , ...
  'YGrid'       , 'off'      , ...
  'XGrid'       , 'off'      , ...
  'GridLineStyle' , '--'       , ...
  'XColor'      , [.05 .05 .05], ...
  'YColor'      , [.05 .05 .05],...
  'LineWidth'   , 2)

set(fh, 'color', 'white');

set(gcf, 'PaperPositionMode', 'auto');