1D Ogata-Banks Solution¶

The figure below depicts the physical system, and the analytical model system at three different times.

Panel A is a depiction of the physical system and the concentration profile along the x-axis at time less than zero. The concentration is zero everywhere.

Panel B is a depection of the physical system and the concentration profile along the x-axis at time equal zero (like the Big Bang!). At x < 0, the concentration is suddenly raised to a value of \(C_0\) everywhere to the left of the origin, x=0. This condition represents a step function input and is a suitable approximation of some upstream source zone that has a constant concentration. The concentration to the right of the origin (x > 0) is still zero.

Panel C is a depiction of the physical system and the concentration profile aling the x-axis at some time greater than zero. The source mass has moved to the right of the origin a distance determined by the species velocity and dispersed along the translational front proportional to the dispersivity in the system.

The analytical solution (Ogata and Banks, 1961) for this situation is:

\[ C(x,t) = \frac{C_0}{2}[erfc(\frac{x-vt}{2\sqrt{D_xt}})+ exp(\frac{xv}{D_x}) \cdot erfc(\frac{x+vt}{2\sqrt{D_xt}})]\]

where

\(C_0\) is the source concentration,
\(D_x = \alpha_L*v+D_d\) is the hydrodynamic dispersion (plus any molecular diffusion),
\(v = \frac{q}{n}\) is the pore velocity

It is also presented in the textbook as equation 6.17

The solution is applicable for porous media flow, where the velocity is the pore velocity (seepage velocity divided by the porosity).

Building a Model¶

This particular model is also all over the internet as an on-line calculator, but sometimes we have need to build our own version. In Jupyter Notebooks running a python kernel we need to script the equation above and structure the equation into a useable function.

Forward define the functions¶

This step is important, the functions must be defined before they are called – in an interpreter, this is usually done at the top of the script. Other scripting languages store the scripts at the end (JavaScript usually keeps scripts at end of the file – it internally promotes then to the top before it runs its JIT bytecode compiler).

In a compiled language, this step is not as necessary (predefinition is, location not so much).

These prototype functions are usually written so that they are organic with respect to their variables, so there is no leakage – in these two functions, the input list is just names, and the output is just a value that must be assigned in the calling script.

#
# prototype ogatabanks function
#
def ogatabanks(c_source,space,time,dispersion,velocity):
    from math import sqrt,erf,erfc,exp  # get special math functions
    term1 = erfc(((space-velocity*time))/(2.0*sqrt(dispersion*time)))
    term2 = exp(velocity*space/dispersion)
    term3 = erfc(((space+velocity*time))/(2.0*sqrt(dispersion*time)))
    print(term3)
    ogatabanks = c_source*0.5*(term1+term2*term3)
    return(ogatabanks)

Concentration Profile¶

The example below uses the function named ogatabanks to compute values for \(x \ge 0\) and \(t \gt 0\), and plot the resulting concentration profile.

#
# example inputs
#
c_source = 1000.0 # source concentration
space = 100.      # how far in X-direction to extend the plot
time  = 50.       # time since release
dispersion = 1.0  # dispersion coefficient
velocity = 1.0    # pore velocity
#
# forward define and initialize vectors for a profile plot
#
how_many_points = 50   
deltax = space/how_many_points
x = [i*deltax for i in range(how_many_points)] # constructor notation
c = [0.0 for i in range(how_many_points)]      # constructor notation
#
# build the profile predictions
#
for i in range(0,how_many_points,1):
    c[i] = ogatabanks(c_source,x[i],time,dispersion,velocity)
#
# Import graphics routines for picture making
#
from matplotlib import pyplot as plt
#
# Build and Render the Plot
#
plt.plot(x,c, color='red', linestyle = 'solid')  # make the plot object
plt.title(" Concentration Profile \n Time: " + repr(time) + " time units \n" + " Dispersion: " + repr(dispersion) + " length^2/time units \n" + " Velocity: " + repr(velocity) + " length/time units \n") # caption the plot object
plt.xlabel(" Distance from Source ") # label x-axis
plt.ylabel(" Concentration        ") # label y-axis
#plt.savefig("ogatabanksplot.png") # optional generates just a plot for embedding into a report
plt.show() # plot to stdio -- has to be last call as it kills prior objects
plt.close('all') # needed when plt.show call not invoked, optional here
#sys.exit() # used to elegant exit for CGI-BIN use

733031437583892e-07
9928852633866985e-07
664089697085712e-08
1435180516621863e-08
631491956652328e-09
9731752900754024e-09
646316074086536e-10
553769516341967e-10
111577818799048e-11
0461915088289236e-11
55962508777167e-12
021255962234891e-13
3618449781240147e-13
9613074980096075e-14
190717543917407e-15
2441921148543639e-15
4038703085471545e-16
464786394576111e-17
971609925696384e-18
368161537187115e-18
2571768119076845e-19
5794976240281e-20
4563071426923016e-21
994442411452461e-22
12585646227532e-22
5239706048321186e-23
982725024512e-24
479331959168028e-25
9798022545929787e-26
442035678959576e-27
821319148997366e-28
0773350070897846e-29
18119084347711e-30
12078247231758e-31
9031145863391725e-32
552964224155404e-33
1082395727792427e-34
61323596624926e-35
1114451161772786e-36
6395123426326041e-37
2234328799099841e-38
77350542614852e-40
046315471890479e-41
004334447363078e-42
5485262910136937e-43
5587073638385598e-44
161241103789536e-46
174351850804583e-47
808453733925338e-48
464826223965468e-49

Concentration History¶

Its also quite reasonable to build a history (concentration versus time at some location) - the script is practically the same, except time is the variable instead of space.

#
# example inputs
#
c_source = 1000.0 # source concentration
space = 50.      # where in X-direction are we
time  = 100.      # how far in T-direction to extend the plot
dispersion = 1.0  # dispersion coefficient
velocity = 1.0    # pore velocity
#
# forward define and initialize vectors for a profile plot
#
how_many_points = 50   
deltat = time/how_many_points
t = [i*deltat for i in range(how_many_points)] # constructor notation
c = [0.0 for i in range(how_many_points)]      # constructor notation

t[0]=1e-5 #cannot have zero time, so use really small value first position in list
#
# build the profile predictions
#
for i in range(0,how_many_points,1):
    c[i] = ogatabanks(c_source,space,t[i],dispersion,velocity)
#
# Import graphics routines for picture making
#
from matplotlib import pyplot as plt
#
# Build and Render the Plot
#
plt.plot(t,c, color='red', linestyle = 'solid')  # make the plot object
plt.title(" Concentration History \n Space: " + repr(space) + " space units \n" + " Dispersion: " + repr(dispersion) + " length^2/time units \n" + " Velocity: " + repr(velocity) + " length/time units \n") # caption the plot object
plt.xlabel(" Time since release ") # label x-axis
plt.ylabel(" Concentration        ") # label y-axis
#plt.savefig("ogatabanksplot.png") # optional generates just a plot for embedding into a report
plt.show() # plot to stdio -- has to be last call as it kills prior objects
plt.close('all') # needed when plt.show call not invoked, optional here
#sys.exit() # used to elegant exit for CGI-BIN use

0
952126631006914e-149
949433113257988e-81
787880039809126e-59
2114989528830613e-47
846411842405321e-41
0404801339498303e-36
1242946645838765e-33
8735664705504999e-31
972385010735985e-30
7941524248040494e-28
9003332723789374e-27
2494563672702272e-26
691873392179496e-26
9432711659535435e-25
267118704691593e-25
1834353814731251e-24
278277449781446e-24
855947703432936e-24
855188239681595e-24
107671698367187e-24
0373234658249127e-23
2395783470276733e-23
395859548981437e-23
49224472995004e-23
5239706048321186e-23
4946615517744525e-23
4141363706883272e-23
2956737967756765e-23
1534510708533014e-23
0005820770300741e-23
479025945460384e-24
034487747459961e-24
724668335212883e-24
577612916653271e-24
6020833685826365e-24
793037343442723e-24
136620622905555e-24
6142490169327249e-24
2056553641217127e-24
909652524634083e-25
519616607890451e-25
7273043524457995e-25
3987035143440443e-25
4242305745001693e-25
7164285941070979e-25
206921284979329e-25
431907953605222e-26
855213715817871e-26
042886119526574e-26

Example 6.1 (Using Near-Field Solution above)¶

From pp. 174-175:

\(K=2.15~\frac{m}{day}\)
\(n=0.1\)
\(\frac{dh}{dx}=0.04\)
\(C_0=1000~ppm\)
\(\alpha_L = 7.5~m\)

Plot a concentration profile at 728 days.

Solution¶

A little preparatory arithmetic is used:

\(v = \frac{1}{n} \cdot K \cdot \frac{dh}{dx} = \frac{1}{0.1} \cdot 2.15 \cdot 0.04 = 0.86~\frac{m}{day}\)
\(D~\approx~\alpha_L \cdot v = 7.5 \cdot 0.86 = 6.45~\frac{m^2}{day}\)

Now just use our program

#
# example inputs
#
c_source = 1000.0 # source concentration
space = 1500.      # how far in X-direction to extend the plot
time  = 728.       # time since release
dispersion = 6.45  # dispersion coefficient
velocity = 0.86   # pore velocity
#
# forward define and initialize vectors for a profile plot
#
how_many_points = 100   
deltax = space/how_many_points
x = [i*deltax for i in range(how_many_points)] # constructor notation
c = [0.0 for i in range(how_many_points)]      # constructor notation
#
# build the profile predictions
#
for i in range(0,how_many_points,1):
    c[i] = ogatabanks(c_source,x[i],time,dispersion,velocity)
#
# Import graphics routines for picture making
#
from matplotlib import pyplot as plt
#
# Build and Render the Plot
#
plt.plot(x,c, color='red', linestyle = 'solid')  # make the plot object
plt.title(" Concentration Profile \n Time: " + repr(time) + " time units \n" + " Dispersion: " + repr(dispersion) + " length^2/time units \n" + " Velocity: " + repr(velocity) + " length/time units \n") # caption the plot object
plt.xlabel(" Distance from Source ") # label x-axis
plt.ylabel(" Concentration        ") # label y-axis
#plt.savefig("ogatabanksplot.png") # optional generates just a plot for embedding into a report
plt.show() # plot to stdio -- has to be last call as it kills prior objects
plt.close('all') # needed when plt.show call not invoked, optional here
#sys.exit() # used to elegant exit for CGI-BIN use

print("Concentration at x = 750, t = 728 :",round(ogatabanks(c_source,750.0,728.0,dispersion,velocity),3)," mg/L")

0432523079828371e-10
70717309820489e-11
2867718381422252e-11
3627157057623855e-12
4447774727577688e-12
673326706246197e-13
4764723869713558e-13
556084612044626e-14
3731536068160933e-14
042046584808178e-15
1620709252130692e-15
2629292401490343e-16
947902056656552e-17
396450351771896e-17
268225936633244e-18
601199364275134e-18
994537631754774e-19
732041328750892e-20
3155483679361915e-20
380375606400359e-21
2208903575350298e-21
705470009804643e-22
854739596363408e-23
2372805546281317e-23
553416215403102e-24
145944007202974e-25
0127350394013818e-25
9463071492538337e-26
65265968833764e-27
694019404311472e-28
1979601371113114e-28
0935012453742212e-29
572542743921763e-30
953234279792017e-31
687191769686667e-32
539254136265097e-32
3882977233686734e-33
6185165517812283e-34
353464766133088e-35
733925828943795e-36
090998528733104e-36
5028164109253448e-37
021360805096524e-38
654828639134684e-39
4047288613574742e-40
263641124293693e-41
21350213321183e-42
224837858710929e-43
257297364206664e-44
261692561449586e-45
183521849997486e-46
967690421098455e-47
0563888692951938e-47
0931910730493809e-48
1046134792236207e-49
0898472911187367e-50
0499317855245331e-51
87633346344873e-53
071283955412949e-54
135407060906072e-55
124042510187317e-56
091294627039051e-57
0854316777857356e-58
145533237797709e-59
2996403490732214e-60
5643997340803997e-61
945974120877961e-62
4418495366806557e-63
0431186534540256e-64
368476948570833e-66
082195554618092e-67
422583350923852e-68
250531637650662e-69
444920869167685e-70
057957853600064e-72
544255950476065e-73
3134655326864715e-74
9335140562981618e-75
1016341565053263e-76
128471066113881e-78
3288259786271926e-79
7654433145548443e-80
141968836714584e-82
622197724400721e-83
281811095300838e-84
0998483621034505e-85
176157413144191e-87
3785006100926542e-88
0671360421094631e-89
674731245539161e-91
9994634247953964e-92
350065142227582e-94
404751978710122e-95
3555013643860988e-96
2690576283758394e-98
9997886885270777e-99
410595383692672e-101
6812641610523584e-102
472034957185388e-104
2671107349634333e-105

9.183521849997486e-46
Concentration at x = 750, t = 728 : 112.838  mg/L

Compare the result to the textbook, which uses the far-field solution for a value of \(100~mg/L\) (how would you modify for the far-field solution?)

A spreadsheet implementation is available below in the references.

Groundwater Transport Phenomena