# DiffusionBias: Bias in Lagrangian Models with Spatially Variable Diffusion (v1.0)

## GETTING STARTED

(c)Ferdi Hellweger, Northeastern University, Boston.

NOTE: The applet requires Java 5 or higher. Java must be enabled in your browser settings. Mac users must have Mac OS X 10.4 or higher. Windows and Linux users may obtain the latest Java from Sun's Java site.

QUICK START: Click setup, then go.

## WHAT IS IT?

[extract from Hellweger and Bucci, 2009]

Microbes can be motile. Chemotaxis by bacteria using flagella and diel vertical migration by phytoplankton using buoyancy regulation are two examples. These active transport behaviors are straightforward to implement in an IBM. However, microbes are also often transported passively by the ambient hydrodynamics. Advection is simulated by simply applying the velocity to the individual’s position, and diffusion can be simulated using a random walk scheme. For a constant diffusivity, the step size is drawn from a normal distribution with standard deviation based on the diffusivity: where z (m) is the position (e.g. depth), R is a random number from a standard normal distribution, K (m2 s-1) is the diffusivity and t (s) is time. However, the diffusivity often varies in space. Mixing in the ocean is driven by wind stress at the water surface and the diffusivity generally decreases with depth. In the metalimnion (thermocline) of a lake, density stratification has a stabilizing effect that leads to lower diffusivity. Application of the simple random walk model to this situation leads to an unrealistic net transport towards areas of lower diffusivity. This can be explained theoretically and numerically, or conceptually by the fact that “over a given time interval, particles at a particular location are influenced by slightly more energetic eddies originating in the area of high diffusion compared with somewhat less energetic eddies from the low diffusion region” (Visser, 1997). The problem can be avoided by transforming the coordinate system so that the diffusivity is constant or by using a different formulation for the step size (Hunter et al., 1993; Visser, 1997; Nagai et al., 2003): where K'(z) is the derivative of K(z). This formulation is different from the simple random walk in that it includes a corrective advection term and the diffusivity is taken at an offset location. Nagai et al. (2003) found that the correction scheme does not work for zero diffusivity and solved this problem by incorporating background diffusivity. Nopass boundaries (e.g. water surface) are typically simulated as reflecting.

## HOW TO USE IT

1. Adjust the slider parameters (see below), or use the default settings.
2. Press the "setup" button.
3. Press the "go" button to begin the simulation.
4. View the bacteria population develop over time.
5. Turn on the "Correction" and observet the effect.

## THINGS TO TRY

Try adjusting the parameters. Use a different pattern for the diffusion. Can the correction remove the bias for all patterns and parameters?

## CREDITS AND REFERENCES

Diffusion Bias References:

Hellweger, F. L., Bucci, V. 2009. A bunch of tiny individuals – Individual-based modeling for microbes (review paper). Ecological Modelling, 220(1):8-22.

Hunter, J.R., Craig, P.D., Phillips, H.E., 1993. On the use of random walk models with spatially variable diffusivity. J. Comput. Phys. 106, 366–376.

Visser, A.W., 1997. Using random walk models to simulate the vertical distribution of particles in a turbulent water column. Mar. Ecol. Prog. Ser. 158, 275–281.

Nagai, T., Yamazaki, H., Kamykowski, D., 2003. A Lagrangian photoresponse model coupled with 2nd-orderclosure. Mar. Ecol. Prog. Ser. 265, 17–30.

NetLogo References:

Wilensky, U. & Reisman, K. (1999). Connected Science: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. International Journal of Complex Systems, M. 234, pp. 1 - 12. (This model is a slightly extended version of the model described in the paper.)

Wilensky, U. & Reisman, K. (in press). Thinking like a Wolf, a Sheep or a Firefly: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. Cognition & Instruction.

## PROCEDURES

```; diffusionbias.nlogo
; ferdi hellweger
; ferdi@coe.neu.edu
; 1/12/2010

globals
[
dt time
dtp ntp
xy xy_u xy_i
xy_p xy_p_u xy_p_i
y_p y_p_u y_p_i y_p_t
xt yt
x2 y2 dx2 dy2 y3 y4
xp
z z_u
a_p_u
v_p_u
maxx minx maxy miny
pcolmax pcolmin pcolmid pcolslope
DmaxPrev DminPrev PatternPrev
D_p_y D_p_t
D_p_u D_p_a D_p_b Dp_p_u
D_b D_b_u Dp_b Dp_b_u D_b3 D_b_u3
nb_u bn bd bn_a bd_a bd_a_p
AnalyticTest mp
]

breed [ bguys bguy ]
patches-own [ D_p Dp_p bn_a_p ]
bguys-own [ ]

to setup

clear-all

set time 0 ;d
set dt 0.01 ;d
set dtp 0.01 ;d
set ntp dt ;d

set xy 10 ;cm
set xy_u ( xy / 100 ) ;m
set xy_i ( max-pxcor - min-pxcor + 1) ;i
set xy_p ( xy / xy_i ) ; cm
set xy_p_u ( xy_p / 100 ) ;m
set xy_p_i 1 ;i
set z 1 ;cm
set z_u ( z / 100 ) ;m
set a_p_u ( xy_p_u * z_u ) ;m2
set v_p_u ( a_p_u * xy_p_u) ;m3
set xp 1e-9
set maxx ( max-pxcor + 0.5 - xp )
set minx ( min-pxcor - 0.5 + xp )
set maxy ( max-pycor + 0.5 - xp )
set miny ( min-pycor - 0.5 + xp )

set pcolmax 17
set pcolmin 12
set pcolmid ( ( pcolmax + pcolmin ) / 2 )
set pcolslope ( ( pcolmax - pcolmin ) / ( Dmax - Dmin ) )

;random-seed 1
set nb_u nb
set-default-shape bguys "circle"
if Initial = "Row"
[
set yt random-pycor
]
if Initial = "Point"
[
set yt random-ycor
set xt random-xcor
]
create-bguys nb_u
[
set color black
set size 0.5
if Initial = "Spread"
[
setxy random-xcor random-ycor
]
if Initial = "Row"
[
setxy random-xcor yt
]
if Initial = "Point"
[
setxy xt yt
]
]

setD

set AnalyticTest false

end

to setD

if not ( Dmax = DmaxPrev ) or not ( Dmin = DminPrev ) or not ( Pattern = PatternPrev )
[
set-current-plot "Diffusion"
set-current-plot-pen "D"
plot-pen-reset
set D_p_y []
set y_p []
set y_p_i n-values ( max-pycor + 1 ) [?]
foreach y_p_i
[
set y_p_t ( ? * xy_p )
if Pattern = "Uniform"
[
set D_p_t ( Dmax - Dmin ) / 2
]
if Pattern = "Linear"
[
set D_p_t Dmin + y_p_t * ( Dmax - Dmin ) / ( xy - xy_p )
]
if Pattern = "Exponential"
[
set D_p_t ( 10 ^ ( ( log Dmin 10 ) + y_p_t * ( ( log Dmax 10 ) - ( log Dmin 10 ) ) / ( xy - xy_p ) ) )
]
if Pattern = "Step"
[
ifelse y_p_t < ( xy / 2 )
[
set D_p_t Dmin
]
[
set D_p_t Dmax
]
]
ask patches with [ pycor = ? ]
[
set D_p D_p_t
]
set D_p_y lput D_p_t D_p_y
set y_p lput y_p_t y_p
plotxy D_p_t y_p_t
]
[
ifelse pycor = max-pycor
[
set D_p_a D_p
]
[
set D_p_a [ D_p ] of patch pxcor ( pycor + 1 )
]
ifelse pycor = min-pycor
[
set D_p_b D_p
]
[
set D_p_b [ D_p ] of patch pxcor ( pycor - 1 )
]
set Dp_p ( D_p_a - D_p_b ) / ( 2 * xy_p )
]

[
set pcolor ( pcolmax - D_p * pcolslope )
]

set DmaxPrev Dmax
set DminPrev Dmin
set PatternPrev Pattern

]

end

to go

setD

[
set D_b [ D_p ] of patch-here
set D_b_u ( D_b * 1.0e-5 / ( 100 ^ 2 ) * 86400 )
set Dp_b [ Dp_p ] of patch-here
set Dp_b_u ( Dp_b * 1.0e-5 / ( 100 ^ 2 ) * 86400 * 100)
set dx2 ( ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u * dt ) )
ifelse Correction
[
set y3 ( ycor + 0.5 * Dp_b_u * dt / xy_p_u )
if y3 > maxy
[
set y3 maxy
]
if y3 < miny
[
set y3 miny
]
set y4 ycor
set ycor y3
set D_b3 [ D_p ] of patch-here
set ycor y4
set D_b_u3 ( D_b3 * 1.0e-5 / ( 100 ^ 2 ) * 86400 )
set dy2 ( Dp_b_u * dt + ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u3 * dt ) )
]
[
set dy2 ( ( random-normal 0.0 1.0 ) * sqrt( 2 * D_b_u * dt ) )
]
set x2 ( xcor + dx2 / xy_p_u )
if x2 > maxx
[
set x2 ( maxx - ( x2 - maxx ) )
]
if x2 < minx
[
set x2 ( minx - x2 )
]
set y2 ( ycor + dy2 / xy_p_u )
if y2 > maxy
[
set y2 ( maxy - ( y2 - maxy ) )
]
if y2 < miny
[
set y2 ( miny - y2 )
]
set xcor x2
set ycor y2
]

if time >= ntp
[

set-current-plot "Bacteria"
set-current-plot-pen "B"
plot-pen-reset
foreach y_p_i
[
set y_p_t ( ? * xy_p )
set bn count bguys-on patches with [ pycor = ? ]
set bd ( bn / xy_p ) ;no/cm
plotxy bd y_p_t
]

if AnalyticTest and Initial = "Row" and Pattern = "Uniform"
[
set mp ( nb_u / xy_i / a_p_u ) ;no/m2
[
set D_p_u ( D_p * 1.0e-5 / ( 100 ^ 2 ) * 86400 )
set y_p_u ( ( pycor - yt ) * xy_p_u )
set bd_a_p ( mp / ( 2 * sqrt ( pi * D_p_u * time ) * exp ( - y_p_u ^ 2 / ( 4 * D_p_u * time ) ) ) ) ;no/m3
set bn_a_p ( bd_a_p * v_p_u )
]
set-current-plot-pen "B_a"
plot-pen-reset
foreach y_p_i
[
set y_p_t ( ? * xy_p )
set bn_a ( sum [ bn_a_p ] of patches with [ pycor = ? ] )
set bd_a ( bn_a / xy_p ) ;no/cm
plotxy bd_a y_p_t
]
]

set ntp ( time + dtp )

]

set time ( time + dt )

end

```