# Lab Transport Modelling
#
rm(list=ls()) # remove all variables from workspace
# Define variables
# Groundwater flow
ne = 0.3 # m^3/m^3 -> effective porosity
gradhx=-0.1 # m/m -> hydraulic gradient in x-direction
gradhy=0 # m/m -> hydraulic gradient in y-direction
Kd=1e-4*60*60*24 # m/d -> hydraulic conductivity
d=10 # m -> aquifer thickness
jwx=-Kd*gradhx / ne # m/d -> flow velocity = filter velocity / porosity
A = 1 * d # m2 -> crossectional area for 1D flow, assuming the aquifer has a width of 1m and d= 10 m
# Transport parameters
Mo=10 # kg/m^3 -> mass added
alpha_L = 1 # m -> longitudinal dispersivity
to=100 # d -> evaluation time
Dxx = alpha_L*jwx # m2/d -> dispersion coefficient in x-direction, calculate from dispersivity and flow velocity
Dyy = 0.1*Dxx # m2/d -> same as above, in y direction, factor as indicated in the problem set
# Estimate location of the peak
xcent=jwx*to # 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
# given that sigma depends on the dispersivity (see lecture notes) and 99.7% of the spread is covered by mu +- 3*sigma,
# I select the plotting area to cover xcent+- 3 sigma
sig3 = 3*sqrt(Dxx*2*to) # m -> only calculate in longitudinal direction
# Near the peak, the solution should be given in the following resolution
xstep=1e-1 # m -> resolution in x direction
ystep=xstep # m -> 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
### Propblem 1a: One dimensional solution
Ct_1d = (Mo/(ne*A)) * (1/(2*sqrt(Dxx*to*pi))) * exp (-(x-jwx*to)^2/(4*Dxx*to))
plot(x,Ct_1d, xlim = c(xcent-sig3,xcent+sig3), xlab="x [m]", ylab= expression(paste("c [kg ",m^{-3},"]")))
### Problem 1b: Two dimensional solution
# Initialize a matrix to receive the solution
Cto_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)) {
Cto_2d[,i]=(Mo/(4*pi*ne*d*to*sqrt(Dxx*Dyy)))*exp((-(x-jwx*to)^2/(4*Dxx*to)) - (y[i]^2/(4*Dyy*to)) )
}
# Plot the solution
# I am setting the plot size automatically to quadratic (3 x 3 inch) in order to achieve equal visual resolution in x and y
# more information on manipulation plots search for help on "par"
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,Cto_2d,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")
###### Problem 1c: Vary alpha_L and alpha_T by an order of magnitude
# initialize a matrix to receive the solution
# decrease D
Dxx_v2 = Dxx/10
Dyy_v2 = Dyy/10
Ct_2d_v2<-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_v2[,i]=(Mo/(4*pi*ne*d*to*sqrt(Dxx_v2*Dyy_v2)))*exp((-(x-jwx*to)^2/(4*Dxx_v2*to)) - (y[i]^2/(4*Dyy_v2*to)) )
}
# Plot the solution
# I am setting the plot size automatically to quadratic (3 x 3 inch) in order to achieve equal visual resolution in x and y
# more information on manipulation plots search for help on "par"
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_v2,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")
###### Problem 1d: Plot breakthroughcurves
# For breakthrough, we consider the concentration distribution at a fixed point over time
xo = 70 # m -> fixed coordinate in x-direction
yo = 0 # m -> fixed coordinate in y-direction
# check in which time range we expect the tracer to arrive. Again, considering that sigma of x
# which is known and now re-arranging for t, we find a quadratic equation f(t) = 0 with the
# following coeficients
a=jwx^2
b=-(18*Dxx+2*jwx*xo)
c=xo^2
t_fast = (-b-sqrt(b^2-4*a*c))/(2*a)
t_slow = (-b+sqrt(b^2-4*a*c))/(2*a)
t_step = 1/24 # h
t=seq(t_fast,t_slow, by=t_step)
# calculate breakthrough curve -> not the scew!
Cxo_1d = (Mo/(ne*A)) * (1/(2*sqrt(Dxx*t*pi))) * exp (-(xo-jwx*t)^2/(4*Dxx*t))
# plot the breaktthough corve
plot(t,Cxo_1d)
