Problem 3a / Headstart

The first part of this script plots the analytical solution. The code is also contained in the R-script “headstart_PBS3.r” which was provided in the last lab session. Note the difference to the analytical solution of Problem Set 1 and 2. There we looked at an instantanous source (a pulse input), here we consider a constant input where the tracer is continously added at the inlet of the pipe (see lecture notes for the difference in the analytical solution and respective plots).

Additional note:

A constant input can be considered as the combination of many instantanous pulses occuring sequentially. Thus, the integral of the pulse response is to be expected for the change of concentration with time and along the flow direction. The response for the pulse input was the Gaussian distribution. Therefore, the response for constant input (many pulses following each other), must be related to the integral of the Gaussian distribution. In fact, it yields an inverse S-shaped function, as shown in the lecture notes, and seems logical. The tracer spreads from zero concentration (far away from the inlet) to the input concentration (near the inlet). In this case, the inflexion point marks the average travel distance since the beginning of the tracer addition, given by

\(x_c = j_{wx}\) * \(t_{out}\).

That “inverse S-shaped” function is decribed with the so-called complementary error function (erfc(x)).

erfc(x) = 1-erf(x)

which depends on the error function (erf(x)). The latter is defined as the integral of a Gauss-type function: erf(x) = \(\frac{1}{\sqrt{\pi}}\) \(\int_{-x}^x\) \(e^{-t^2}\) dx .

But now, let’s get started:

# This will be a one dimensional model. Before you start, please install an additional library: "NORMT3" Do this by chosing the tab "packages" in the bottom right window, click "install" and type NORMT3 in the "Packages" space.

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
Cbound = 10           # kg/m^3 -> concentration boundary at the left side
alpha_L = 1e-3        # m -> longitudinal dispersivity
to=5                  # d -> evaluation time
Dxx = alpha_L*jwx     # m2/d -> dispersion coefficient in x-direction, calculate from dispersivity and flow velocity


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

# 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 

# The solution should be given in the following resolution 
xstep=0.01         # m -> resolution in x direction
ystep=xstep        # m -> same resolution in y direction

# find the location where the concentration front is located 
x=seq(xcent-sig3,xcent+sig3, by=xstep)  

# One-dimensional solution for constant input. In order to be able to use the error function, we need to install a new package NORMT3: 
library("NORMT3")
## NORMT3: Evaluates erf, erfc, Faddeeva functions and Gaussian/T sum densities
## Copyright: Guy Nason 2005-2012
Cana=(Cbound/2)*erfc((x-jwx*to)/(2*sqrt(Dxx*to)));

plot(x,Cana, xlim = c(xcent-sig3,xcent+sig3), xlab="x [m]",  ylab= expression(paste("c  [kg ",m^{-3},"]")))

Numerical solution

The Courant number is given (Co = 0.5), and this defines the time step.

# preparation: Calculate the time step such that the Courant number < 1
CourantNo=0.5
deltax=xstep
#deltat=0.5*deltax^2/(jwx*alpha_L)
deltat=CourantNo*deltax/jwx # days
A side remark on numerical errors:

We can calculate the grid Peclet number, and we will realize: It is much too small. Remember, in order to avoid numerical dispersion, Peg should be > 1/2.

(Peg=Dxx*deltat/deltax^2)  
## [1] 0.05

Thus, it is no suprise that we will run into numerical dispersion (see plot below). On the other hand, increasing Peg will also increase Co, which however will enhence the risk to run into oscillations. In this example, Co may be elevated to around Co=0.8 without yielding oscillations, but this will not yet avoid numerical dispersion completely. There is no perfect solution.

Here is the explicit solution:

# initialize arrays
tDim=round(to/deltat)
xDim=round(max(x)/deltax) - 1

DT=Dxx*deltat/(deltax^2) # Factor for dispersion
AT=jwx*deltat/deltax     # Factor for advection

Cexp = matrix(0, nrow = xDim, ncol = tDim+1) # each x one row, each t onw col
Cexp[,1]=Cbound

# external loop: loop through time
for (t in 2:tDim){  # n, cols, time
# boundary condition at the left boundary (C(x=0)=C_bound)
j=1
Cexp[j,t+1]=Cbound

# internal loop: through space
for (j in 2:(xDim-1)){ # j, rows, space
Cexp[j,t+1]=Cexp[j,t]+DT*(Cexp[j-1,t]-2*Cexp[j,t]+Cexp[j+1,t])-AT*(Cexp[j,t]-Cexp[j-1,t])
}

}

## compare to the analytical solution:
plot(x,Cana, xlim = c(xcent-sig3,xcent+sig3), xlab="x [m]",  ylab= expression(paste("c  [kg ",m^{-3},"]"))) # analytical solution
lines(x, Cexp[(xDim-length(x)+1):xDim,tDim],type="l") # numerical solution