Prologue

# +++ Solution to Problem Set 1 Transport modelling (analytical solution) +++
# Author: Anke Hildebrandt
# Date: 22.05.2018

Definition of Variables

Note: The flow velocity is found from the Darcy velocity: \[ j_{w,F} = K \frac{d h} {d x} \] and the actual flow velocity is: \[ j_{wa,x} = \frac{j_{w,F}}{n_e}\]

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  
jwx=-Kd*gradhx / ne     # m/d -> flow velocity = filter velocity / porosity
d=10                    # m -> aquifer thickness
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*60*60*24               # d -> evaluation time
(Dxx = alpha_L*jwx)      # m2/d -> dispersion coefficient in x-direction, calculate from dispersivity and flow velocity
## [1] 3.333333e-05
Dyy = 0.1*Dxx          # m2/d -> same as above, in y direction, factor as indicated in the problem set

Problem 1a - One dimensional solution

Estimate the location and width of the peak: The centre of the mass distribution is given by the flow velocity and the time elapsed.

The centre of mass is than calculated as \(x_{cent} = j_{wa,x} * t_1\)

# Estimate location of the peak
(xcent=jwx*to)      # m -> centre of the distribution in x direction, calculate from advection
## [1] 288

The width of the peak is determined by the dispersion, with $ _x = $ Given that the spread around the centre depends on the dispersivity (see lecture notes) and 99.7% of the spread is covered in a normal distribution by \(\mu \pm 3 \cdot \sigma\), I select the plotting interval to \(x_{cent} \pm 3 \cdot \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
x=seq(xcent-sig3,xcent+sig3, by=xstep)  # these are the x-points the solution will be calculated on

Now I am ready to calculate and plot the 1-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 analytical solution

We use the equation given in the lecture (instantanous input) for two dimensions. For this, I need to define also the y-axis, and I give it the same resolution as for x. Also, I am setting the expected centre of the mass distribution in y-direction. It is located at 0 m, since the ground water flow is only in x-direction and therefore also advection occurs only in x-direction. With the injection taking place at x=0,

ycent = 0           # m -> centre of the distribution in y direction 
ystep=xstep

y=seq(ycent-sig3,ycent+sig3, by=ystep)  # these are the x-points the solution will be calculated on

Now I am ready to calculate the two dimenstional solution. I calculate it in slices, that it for each increment of y the mass at all x. And I store the calculated values sequentially in a mstrix. That matrix needs to be intialized first (we have seen that already in R-intro)

# 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*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 -> for more information on manipulation of 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,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")

Problem 1c - Vary dispersivities

I am hypothesize that the spread of the centre of mass should be located at the same spot - since advection determines the mean location of the peak. I also expect the spread to be smaller, since the variance of the analytical solution depends on dispersivity.

# 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

Again, I am setting the plot to be quadratic, as above.

opar <- par()   # I am saving the current parameters as a backup
par(pin=c(3,3)) # setting the plot size
contour(x,y,Ct_2d_v2,xlim = c(xcent-sig3,xcent+sig3), ylim=c(ycent-sig3,ycent+sig3), xlab="x [m]", ylab="y [m]")