Prologue

rm(list=ls()) # remove all variables from workspace

# Define variables

# Groundwater flow
ne = 0.3                # m^3/m^3
gradhx=-0.1             # m/m
gradhy=0                # m/m
Kd=1e-4*60*60*24        # m/d 
d=10                    # m
jwx=-Kd*gradhx / ne     # m/d -> Darcy * porosity
A = 1 * d              # m2 -> for 1D flow, assuming the aquifer has a width of 1m and d= 10 m

# Transport parameters
Mo=10                  # kg
alpha_L = 1            # m
tout=100               # d
Dxx = alpha_L*jwx      # m2/d -> calculate from dispersivity and flow velocity
Dyy = 0.1*Dxx          # m2/d -> as indicated


# Estimate location of the peak
xcent=jwx*tout      # m -> centre of the distribution in x direction, calculate from advection
ycent = 0           # m -> centre of the distribution in y direction

Estimate the width of the distribution: Remember that the spread of the distribution (\(\sigma\)) depends on the dispersivity and travelling time (see lecture notes). Also, we know that 99.7% of the spread is covered by \(\mu\) +- 3* \(\sigma\), where \(\mu\)= x_cent= \(j_{wx} \cdot t_{out}\). Thus, select the plotting area to cover x_cent +- 3 * \(\sigma\)

sig3 = 3*sqrt(Dxx*2*tout) # see lecture note, only calculate in longitudinal direction

# The solution should be given in the following resolution 
xstep=1e-1         # resolution in x direction
ystep=xstep        # same resolution in y direction

# select the area comprising 99.7 % of the concentration
x=seq(xcent-sig3,xcent+sig3, by=xstep)  # these are the x-points the solution will be calculated on
y=seq(ycent-sig3,ycent+sig3, by=ystep)  # these are the y-points the solution will be calculated on

Two dimensional analytical solution

I am calculating again the 2-dimensional analytical solution, just in order to be able to compare it to the particle tracking below.

# initialize a matrix to receive the solution
Ct_2d<-matrix(NA,length(x),length(y)) # length(x) = number of rows, length(y) = number of columns

# calculate the solution
for(i in 1:length(y)) {
  Ct_2d[,i]=(Mo/(4*pi*ne*d*tout*sqrt(Dxx*Dyy)))*exp((-(x-jwx*tout)^2/(4*Dxx*tout)) - (y[i]^2/(4*Dyy*tout)) )
  }

opar <- par()   # I am saving the current parameters as a backup
par(pin=c(3,3)) # now I am changing the plot parameters on the plot size
# the following is the plot command, it will use the paramers I set above and therefore the plot will be quadratic
contour(x,y,Ct_2d,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")

Patricle tracking

Use the implementation shown in the lecture.

I am stepping each particle one day at a time. Each step consists of a deterministic (advection) and random (dispersion) component. Steps in x and y direction are calculated separatly.

np = 2.5e+4 # number of particles to be implemented in the field 
tstep = 1 # days 
ntime = tout/tstep # number of time steps to be modelled

# initial conditions for particles (all start at x=0, y=0)
pos_x_p = rep(0,np) # this is an array that contains np columns
pos_y_p = rep(0,np) # this is an array that contains np columns

for (ti in 1:ntime) {
  delta_x_p = jwx*tstep+sqrt(Dxx*2*tstep)*rnorm(np,0,1) # step all np particles in x direction
  delta_y_p = 0 + sqrt(Dyy*2*tstep)*rnorm(np,0,1) # step all np particles in y direction
  pos_x_p = pos_x_p + delta_x_p # new x position for all np particles
  pos_y_p = pos_y_p + delta_y_p # new y position for all np particles
}

# plot particle tracking solution at the end
par(pin=c(3,3))
plot(pos_x_p, pos_y_p, cex=0.05, xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")

Optional: Compare particle tracking and analytical solution

Test whether the particle tracking yields the same result as the analytical solution. One option to do this is to calculate the mass distribution obtained from the analytical solution and particle tracking and compare, whether they give the same result. For matters of simplicity, I will not compare the 2D solution directly, but the mass distribution along the x and y axis.

# Initialize a vector which will contain the number of particles contained in each window (or section) of size xstep along the x-axis.
counts_1D_x = matrix(data=0,nrow=length(x), ncol=1)

# count the particles falling in a certain window along x direction
for (i in 1:length(x)) {
  counts_1D_x[i] = length(which(abs(x[i]-pos_x_p)<xstep/2)) 
}
 
# To convert counts to concentration we have to calculate, which mass corresponds to each particle. We added np patricles to simulate Mo kg of added tracer, thus each particle corresponds to 
m_p = Mo/np # kg

# The 1D mass distribution of the average concentration along the x-axis estimated from particle tracking is
m_pt_1d_x = counts_1D_x*m_p  

# The  corresponding mass distribution based on the analytical solution is
m_anal_1d_x = rowSums(Ct_2d*xstep*ystep*d*ne) 
# note 1: concentration * volume = mass
# note 2: rowsum means a sum along the y axis and yields the mass distribution along the x axis

# Overlay the two graphs to see that they do indeed overlap pretty well!
par(pin=c(3,3))
plot(x,m_anal_1d_x, col="red", ylab="Mass [kg]", xlab="x [m]", main="along x-axis")
points(x,m_pt_1d_x, cex=0.05)

# the same can be done in the y-direction:
counts_1D_y = matrix(data=0,nrow=length(y), ncol=1)

# count particles falling in a certain window along y direction
for (i in 1:length(y)) {
  counts_1D_y[i] = length(which(abs(y[i]-pos_y_p)<ystep/2)) 
}

m_pt_1d_y = counts_1D_y*m_p  #  mass distribution rom particle traking, y-direction

# corresponding mass distribution from the analytical solution
m_anal_1d_y = colSums(Ct_2d*xstep*ystep*d*ne) 

plot(m_anal_1d_y,y, col="green", xlab="Mass in [kg]", ylab="y [m]",main="along y-axis")
points(m_pt_1d_y,y, cex=0.5)

Combined plot

Using a function called “layout” (use help or google to learn more), the plots can be arranged in an intuitive fashion.

layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE)) # this arranges the figures
contour(x,y,Ct_2d,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]") 
points(pos_x_p, pos_y_p, cex=0.02)

plot(m_anal_1d_y,y, col="green", xlab="Mass in [kg]", ylab="y [m]")
points(m_pt_1d_y,y, cex=0.5)

plot(x,m_anal_1d_x, col="red", ylab="Mass in [kg]", xlab = "x [m]")
points(x,m_pt_1d_x, cex=0.5)