Modeling Advection by Particle Tracking

Modeling Advection by Particle Tracking#

Later on we will use USGS-MOC (which incorporates particle tracking to define characteristic trajectories), so we will illustrate particle tracking keeping MOC in mind.

  • Particle or front tracking is typically performed using special software.

  • It can be performed using a spreadsheet.

  • The spreadsheet exercise is useful to illustrate the principles involved in particle tracking calculations.

  • Particle tracking with reactions is very computationally intensive and is beyond practical application in a spreadsheet; such computations are used to design space shuttles, helicopters and other things where the knowledge gained justifies the costs

The first step is to compute the velocity field (if it is steady even better, if not then every time increment it is recomputed)

  • If analytical functions are available for the velocity field then tracking is relatively easy.

  • Usually the velocity field is determined numerically at discrete points in space, and this is the situation of interest.

  • The interpolation schemes in common use are simple; simple, simple-linear, and multi-linear schemes.

  • Only the simple-linear scheme preserves cell-by-cell mass balances.

The figure below illustrates simple linear interpolation.

  • Typical computational grid for heads.

  • Arrows are the interfacial fluxes.

  • The simple scheme assigns the top and left flux to \((x_1,y_1)\).

  • The simple scheme assgins the right and bottom flux to \((x_2,y_2)\).

The particle velocity is determined by position of the particle relative to the velocity grid.

Uses the same grid as the head scheme. Velocity is the distance weighted average of the cell that the particle occupies.

From the above figure, the particle’s velocities (x- and y- directions) are:

  • \(u_p = \frac{1-\delta x}{\Delta x}u(x_1,y_1)+\frac{\delta x}{\Delta x}u(x_2,y_2)\)

  • \(v_p = \frac{1-\delta y}{\Delta y}v(x_1,y_1)+\frac{\delta y}{\Delta y}v(x_2,y_2)\)

Multi-Linear Interpolation#

Higher order schemes produce smoother velocity fields at the expense of cell mass balances and computational ease. The USGS-MOC model uses a bi-linear scheme where the velocities at the four corners of the occupied cell are used. When transient flow fields occur, averaging in time is also used. The differences in the schemes are hard to detect when the grid spacing is small and the flow field is smoothly varying.

Spreadsheet Approach#

To illustrate particle tracking the simple velocity scheme is used.
Extension to higher order schemes is straight forward.

Illustrate with simple scheme. Consider:

  • Large rectangles represent the velocity grid.

  • Circles represent the geometric location where velocity is known.

  • Small rectangle represents the particle that we wish to track.

Cell Indexing

  • Each cell represents a grid location in the velocity field. Thus each cell has a unique row and column index.

  • Each cell centroid also has a unique geomteric (x,y) location.

  • The particle in the figure is located in cell named: Col_1,Row_2.

  • The cell is located at position: (X_1,Y_2).

  • The particle position is (XP,YP).

Locating the Particle

  • At the start of a time-step:

  • particle position is known.

  • cell positions are known.

  • cell that the particle occupies is unknown.

  • Construct a distance table

  • The distance from each cell to the particle is calculated and stored in a table.

  • Search the table,find the cell nearest the particle.

  • The cell coordinates of the smallest distance in the table is determined

Locating in Excel

  • The spreadsheet function that finds the value in an array (rectangular area of cells), given the position in the array to search is the function INDEX(array,row_index,column_index)

  • The spreadsheet function that can find the position in an array where a particular value appears is the function MATCH(value,array,type)

INDEX Function

  • INDEX(array,row_index,column_index) array is the location of the rectangular area of cells to search (eg. A3:C6). row_index is number of rows down from the starting row to search. column_index is the number of columns across from the starting column to search.

MATCH Function

  • MATCH(value,array,type) value is the numerical value to search for in the array. array is the location of the rectangular area of cells to search (eg. A3:C6). type is the type of match to use. type=0 means exact matching.

Using the Functions

  • The INDEX function allows us to select the correct values of velocity if we know which cell the particle resides in.

  • The MATCH function allows us to compare values in an array and determine the position in the array that these values are found.
    Thus the MATCH function lets us search a distance table, find the cell center nearest the particle, and then use the index to find the correct velocity.

Moving the Particle

Once the cell containing the particle is identified, the particle is assigned the velocity values for that cell. The particle is then “moved” by the simple kinematic calculation:

  • \(x_p(t+\Delta t)=x_p(t)+u_p(t)\Delta t\)

  • \(y_p(t+\Delta t)=y_p(t)+v_p(t)\Delta t\)

Illustrative Example

An example with 3 particles is shown below - the spreadsheet is kind of big because iterative computations are intentionally avoided.

Closer view for particle 1:

Closer view for particle 2:

Closer view for particle 3:

Summary

  • Particle tracking is a tool to determine the position of a fluid particle in a flow field.

  • A two-step approach is required:

  • Determine particle velocity

    • Locate the particle relative to known velocity locations.

    • Assign the velocity to the particle based on an interpolation scheme.

  • Move the particle.

  • All particle tracking programs use some sort of this two-step logic.

  • Extension to a 3rd spatial dimension is straightforward; boundaries are tricky - ray tracing and folding back across a boundary is a simple approach to manage a need to keep particles in a computational domain; alternatively

Spreadsheet Examples#

  1. Cleveland, T.G. (1997) Particle Tracking Spreadsheet (10X10) (Excel)

  2. Cleveland, T.G. (1997) Particle Tracking Spreadsheet (20X20) (Excel)

  3. Cleveland, T.G. (1997) Particle Tracking Spreadsheet (5X5 same as example above) (Excel)

Performing the calculations in a python script is illustrated in the next section.

Note

You will see that the python script is simpler to construct than the spreadsheet (in terms of weird functions), but takes a shit-ton more code!

Example of Particle Tracking#

Here we illustrate using particle tracking to address a homework problem.

The figure below shows a piezometric map for a shallow sand aquifer. The hydraulic conductivtiy is estimated to be \(1.5 \times 10^{-2}~\frac{cm}{s}\), the saturated thickness is 40 feet, and the effective porosity is 0.3.

Determine:

  1. Which well is expected to be the most contaminated.

  2. The groundwater velocity and seepage velocity across the plume.

  3. The duration that the source has been contaminating the aquifer (neglect dispersion, diffusiom, and adsorption).

  4. The flow rate across the plume.

  5. An explaination for contamination upgradient of the source zone.

Step 1. Need to Digitize the Head Map#

Use G3DATA or something similar to digitize the head map.

Using G3DATA we will need a coordinate system, so thats superimposed on

Next you need to open an instance of G3DATA. If using Windows the system should give you a windows manager, if using linux you need an X-window desktop.

Open G3DATA#

Observe that I get the RandR warning (Resize and rotate) - this is assocated with my VNC server window, and the software still is working. Given I am running on a Raspberry Pi over the internet, ill live with the warning.

Also notice in the program call the -scale 0.1 this is a trial-and-error to rescale the image so I can see what I am doing, so open and close a few times until you get the view window into a workable size.

Set View == Zoom#

Next we will need to set the view window to Zoom area. This is needed to write the level set files that we are building.

Set Filename#

Here we set the file name for each level set. We will digitize each contour line and save its coordinates to a file, then we conactenate these files after adding the Z coordinate into a single file for gridding and plotting.

Set XY Coordinates#

Here we must set the X and Y coordinates of the lower left,lower right (X) and lower left upper left (Y) to establish how we will be measuring things.

Digitize the Contours#

Now we digitize a contour line

when done, we write the points to a file. Here is what is in the tile.

Then clear the points (in the G3DATA program), change the file name and digitize the next line. Repeat until all the lines are complete.

Then open each file and add the Z-coordinate. Lastly join each file into a single file for further processing.

Here is what my version looks like, the file is Fig5.18-LevelSets.png.dat its an ASCII file that I can read and make a contour plot, but more importantly I can grid the data to estimate x- and y- components of velocity from the given value of hydraulic conductivity.

X-Easting	Y-Northing	Z-Elevation
-0.365791375238  0.00426936718372 5.5
7.44799110264  34.0524726574 5.5
30.6624750375  65.7872126049 5.5
...
many lines
...
353.704129652  725.957058705 5.5
186.509893725  0.312464310759 5.75
195.371098732  37.1933927273 5.75
202.181166327  65.9190293167 5.75
206.765265923  96.4487403766 5.75
...
many lines
...
486.889653936  722.268859129 5.75
362.749010841  -18.4586090188 6.0
367.376818083  9.58126046662 6.0
...
many lines
...
595.939191403  718.506746134 6.0
502.803399911  -17.6039350757 6.25
512.271308678  26.3822876764 6.25
520.432150683  61.4932965599 6.25
525.047470026  90.2445493523 6.25
...
many lines
...
708.769440297  624.373336948 6.25
637.755441575  -17.7565649526 6.50
649.674100524  11.6206838032 6.50
...
many lines
...
885.881149356  618.394355043 6.50
901.700715759  -11.9467567219 6.75
909.60191353  17.1217634194 6.75
...
many lines
...
976.769639356  628.357723872 6.75

Step 2 Check the Digitization Results#

Here we will now produce a contour map from the digitized map, adjust sizes and overlay on the original map to get an idea of how good we did.

Also will use trial-and-error to find a useable hull within the contour plot where we can capture heads and use for a particle tracking effort - we have to find a rectangle with non-null interpolation values; easiest way (my opinion) is to try indices until can find a reasonable rectangle then sample from that rectangle. The process is documented below.

Note

The interpolation algorithm does not estimate vales outside the convex hull defined by the observation data. So the box below is by trial-and-error to stay within the convex hull; the interpolator in this region produces useable estimates for computing the heads and velocities from just the mapped information>

# CCMR from ENGR-1330:
# http://54.243.252.9/engr-1330-webroot/8-Labs/Lab07/Lab07.html
# https://clouds.eos.ubc.ca/~phil/docs/problem_solving/06-Plotting-with-Matplotlib/06.14-Contour-Plots.html
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
# https://stackoverflow.com/questions/332289/how-do-you-change-the-size-of-figures-drawn-with-matplotlib
# https://stackoverflow.com/questions/18730044/converting-two-lists-into-a-matrix
# https://stackoverflow.com/questions/3242382/interpolation-over-an-irregular-grid
# https://stackoverflow.com/questions/33919875/interpolate-irregular-3d-data-from-a-xyz-file-to-a-regular-grid
import pandas
my_xyz = pandas.read_csv('Fig5.18-LevelSets.png.dat',sep='\t') # read an ascii file already prepared, delimiter is tabs
#my_xyz = pandas.read_csv('XYZSomewhereUSA.txt',sep=' ') # read an ascii file already prepared, delimiter is tabs
my_xyz = pandas.DataFrame(my_xyz) # convert into a data frame
#print(my_xyz) #examine the dataframe
import numpy 
import matplotlib.pyplot
from scipy.interpolate import griddata
# extract lists from the dataframe
coord_x = my_xyz['X-Easting'].values.tolist()
coord_y = my_xyz['Y-Northing'].values.tolist()
coord_z = my_xyz['Z-Elevation'].values.tolist()
coord_xy = numpy.column_stack((coord_x, coord_y))
# Set plotting range in original data units
lon = numpy.linspace(min(coord_x), max(coord_x), 100)
lat = numpy.linspace(min(coord_y), max(coord_y), 100)
X, Y = numpy.meshgrid(lon, lat)
# Grid the data; use cubic spline interpolation (other choices are nearest and linear)
Z = griddata(numpy.array(coord_xy), numpy.array(coord_z), (X, Y), method='cubic', fill_value = 'nan')
# Build the map
print("Indices of black box on overlay below:")
ixl = 9
ixh = 90
iyl = 27
iyh = 62
print("xyz lower left corner ",X[iyl][ixl],Y[iyl][ixl],Z[iyl][ixl])
print("xyz lower right corner",X[iyl][ixh],Y[iyl][ixh],Z[iyl][ixh])
print("xyz upper left corner ",X[iyh][ixl],Y[iyh][ixl],Z[iyh][ixl])
print("xyz upper right corner ",X[iyh][ixh],Y[iyh][ixh],Z[iyh][ixh])
xxl = X[iyl][ixl]
xxh = X[iyh][ixh]
yyl = Y[iyl][ixl]
yyh = Y[iyh][ixh]
flag=True
if flag:
    matplotlib.pyplot.rcParams["figure.figsize"] = [10.00, 8.00]
    matplotlib.pyplot.rcParams["figure.autolayout"] = True
    im = matplotlib.pyplot.imread("Fig5.18-Coordinates.png") # base image
fig, ax = matplotlib.pyplot.subplots()
if flag: 
    im = ax.imshow(im, extent=[-25, 1000, -125, 800])# sets X and Y plot window of basemap
#fig.set_size_inches(14, 7)
levels=[5.50,5.75,6.00,6.25,6.50,6.75]
CS = ax.contour(X, Y, Z, levels, linewidths=3)
ax.clabel(CS, inline=2, fontsize=12)
ax.set_title('Contour Plot from Gridded Data File')
ax.set_xlim([0,1000])
ax.set_ylim([0,800])
ax.plot([xxl,xxh],[yyl,yyl],color="black")
ax.plot([xxl,xxh],[yyh,yyh],color="black")
ax.plot([xxl,xxl],[yyl,yyh],color="black")
ax.plot([xxh,xxh],[yyl,yyh],color="black");

Indices of black box on overlay below:
xyz lower left corner  90.9090909090909 200.0 5.478952367771485
xyz lower right corner 909.090909090909 200.0 6.7224489855441645
xyz upper left corner  90.9090909090909 588.8888888888889 5.302942021834025
xyz upper right corner  909.090909090909 588.8888888888889 6.663235234576626
../../_images/1a9f6dda9e0b4aa1354015bd7bbc4da23c85c4d26106035c451571e0528dbbda.png 
# === Load Required Libraries ===
import pandas
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri

# === Load Data from ASCII File ===
# This file is expected to have tab-separated columns: X-Easting, Y-Northing, Z-Elevation
my_xyz = pandas.read_csv('Fig5.18-LevelSets.png.dat', sep='\t')
my_xyz = pandas.DataFrame(my_xyz)

# === Extract Columns into Arrays ===
coord_x = my_xyz['X-Easting'].values
coord_y = my_xyz['Y-Northing'].values
coord_z = my_xyz['Z-Elevation'].values

# === Create Triangulation from Scattered Points ===
triang = tri.Triangulation(coord_x, coord_y)

# === Define Bounding Box Indices in Grid Coordinates (manually set) ===
# These are used to extract corners for annotation
ixl = 9   # x-index low
ixh = 90  # x-index high
iyl = 27  # y-index low
iyh = 62  # y-index high

# === Compute Bounding Box Corners for Annotations ===
# Because there is no gridding here, we must derive directly from coordinate logic
# We'll define corners directly using data range rather than relying on mesh indices
# Here we compute approximate bounding box from coordinate values
xxl, xxh = 100, 800  # manual adjustment as needed
yyl, yyh = 200, 600

# === Prepare Plotting Environment ===
fig, ax = plt.subplots(figsize=(10, 8))

# === Overlay Background Image (optional) ===
im = plt.imread("Fig5.18-Coordinates.png")  # base map image
ax.imshow(im, extent=[-25, 1000, -125, 800])  # extent must match image bounds

# === Plot Contour Lines from Triangulated Data ===
contour_levels = [5.50, 5.75, 6.00, 6.25, 6.50, 6.75]
CS = ax.tricontour(triang, coord_z, levels=contour_levels, linewidths=2)
ax.clabel(CS, inline=2, fontsize=12)

# === Annotate Bounding Box (Black Box Over a Region of Interest) ===
# Draw rectangle from (xxl, yyl) to (xxh, yyh)
ax.plot([xxl, xxh], [yyl, yyl], color="black")  # bottom edge
ax.plot([xxl, xxh], [yyh, yyh], color="black")  # top edge
ax.plot([xxl, xxl], [yyl, yyh], color="black")  # left edge
ax.plot([xxh, xxh], [yyl, yyh], color="black")  # right edge

# === Customize Plot Appearance ===
ax.set_title('Triangulated Contour Plot with Overlay and Bounding Box')
ax.set_xlabel("X-Easting")
ax.set_ylabel("Y-Northing")
ax.set_xlim([0, 1000])
ax.set_ylim([0, 800])
ax.set_aspect("equal")  # Maintain 1:1 aspect ratio
ax.grid(True)

# === Show the Final Plot ===
plt.show()
../../_images/3531376059ce7056d055c21093a3f8660287564b77fb6c7271965eead507fa09.png

Contour Plot Using Triangulated Data#

This example demonstrates how to:

  • Generate a contour plot from irregularly spaced data using triangulation

  • Overlay a bounding box to highlight a region of interest

  • Use an interpolation tool to estimate elevation (Z) at an arbitrary location

  • (Optional) Display a base map beneath the contours

Later in this course, you will use the interpolation tool in a particle tracking exercise to simulate flow along hydraulic gradients.

The dataset Fig5.18-LevelSets.png.dat should be located in your working directory and contains X, Y, and Z columns for Easting, Northing, and Elevation.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from matplotlib.tri import LinearTriInterpolator

# === Load Data ===
my_xyz = pd.read_csv('Fig5.18-LevelSets.png.dat', sep='\t')
coord_x = my_xyz['X-Easting'].values
coord_y = my_xyz['Y-Northing'].values
coord_z = my_xyz['Z-Elevation'].values

# === Triangulate for Irregular Grid ===
triang = tri.Triangulation(coord_x, coord_y)

# === Optional Base Map Image ===
use_basemap = True
if use_basemap:
    plt.rcParams["figure.figsize"] = [10.00, 8.00]
    plt.rcParams["figure.autolayout"] = True
    im = plt.imread("Fig5.18-Coordinates.png")

fig, ax = plt.subplots()

if use_basemap:
    ax.imshow(im, extent=[-25, 1000, -125, 800])  # Adjust extent to match coordinates

# === Contour Plot via Triangulation ===
levels = [5.50, 5.75, 6.00, 6.25, 6.50, 6.75]
CS = ax.tricontour(triang, coord_z, levels=levels, linewidths=2)
ax.clabel(CS, inline=True, fontsize=12)

# === Bounding Box Overlay ===
ixl, ixh = 9, 90
iyl, iyh = 27, 62
lon = np.linspace(min(coord_x), max(coord_x), 100)
lat = np.linspace(min(coord_y), max(coord_y), 100)
X, Y = np.meshgrid(lon, lat)

xxl = X[iyl][ixl]
xxh = X[iyh][ixh]
yyl = Y[iyl][ixl]
yyh = Y[iyh][ixh]

ax.plot([xxl, xxh], [yyl, yyl], color="black")
ax.plot([xxl, xxh], [yyh, yyh], color="black")
ax.plot([xxl, xxl], [yyl, yyh], color="black")
ax.plot([xxh, xxh], [yyl, yyh], color="black")

# === Interpolator Block (for future particle tracking logic) ===
interpolate_flag = True
if interpolate_flag:
    interpolator = LinearTriInterpolator(triang, coord_z)

    # Define a test point
    x_query = 400
    y_query = 200
    z_est = interpolator(x_query, y_query)
    print(f"Interpolated Z at ({x_query:.2f}, {y_query:.2f}) = {z_est:.3f}")

    ax.plot(x_query, y_query, 'ro')  # Red dot
    ax.annotate(f"{z_est:.2f}", (x_query + 10, y_query + 10), color='red', fontsize=10)

# === Final Plot Settings ===
ax.set_title("Triangulated Contour Plot with Bounding Box and Interpolation Point")
ax.set_xlim([0, 1000])
ax.set_ylim([0, 800])
ax.set_xlabel("X - Easting")
ax.set_ylabel("Y - Northing")
plt.grid(True)
plt.show()

Interpolated Z at (400.00, 200.00) = 6.001
../../_images/1aab87fedb74098fcfb98a879631c52186924027be3763ada799680d958e5703.png

This is acceptable. Now we will extract velocities from the gridded data files - using the “box”. Recall to use the spreadsheet version of the particle tracker we need the \(U\) and \(V\) fields from the head distribution.

We will use the supplied hydraulic conductivity and porosity

  • \(K~\approx~ 1.5 \times 10^{-2}~\frac{cm}{s} \cdot \frac{1~in}{2.54~cm} \cdot \frac{1~ft}{12~in} \cdot \frac{86400~s}{1~day} = 42.52~\frac{ft}{day}\)

  • \(n~\approx~ 0.30 \)

And apply Darcy’s law as

  • \(\textbf{u(x,y)} = -\frac{K}{n}\cdot \frac{\Delta h}{\Delta x}\)

  • \(\textbf{v(x,y)} = -\frac{K}{n}\cdot \frac{\Delta h}{\Delta y}\)

Conductivity = 42.52 # feet/day
Porosity = 0.30 # voids/bulk
nrows = iyh - iyl
ncols = ixh - ixl
#Zero vectors for the velocity fields
U = [[0 for j in range(ncols)] for i in range(nrows)]
V = [[0 for j in range(ncols)] for i in range(nrows)]
#Zero vectors for realinged mesh grid
XG = [[0 for j in range(ncols)] for i in range(nrows)]
YG = [[0 for j in range(ncols)] for i in range(nrows)]
for j in range(ncols):
    for i in range(nrows): #range(nrows)
        U[i][j] = -1.0*(Conductivity/Porosity)*(Z[iyl+i][ixl+j]-Z[iyl+i][ixl+j-1])/(X[iyl+i][ixl+j]-X[iyl+i][ixl+j-1]) #- K dh/dx
        V[i][j] = -1.0*(Conductivity/Porosity)*(Z[iyl+i][ixl+j]-Z[iyl+i-1][ixl+j])/(Y[iyl+i][ixl+j]-Y[iyl+i-1][ixl+j]) #- K dh/dy
        XG[i][j] = X[iyl+i][ixl+j]
        YG[i][j] = Y[iyl+i][ixl+j]

Suppose we want a plot of the vector field just to check our work, a nice tool is quiver and/or streamplot which are demonstrated below

# convert to numpy arrays
UU = numpy.asarray(U)
VV = numpy.asarray(V)
XX = numpy.asarray(XG)
YY = numpy.asarray(YG)

UU[30][30]

np.float64(-0.2250452891747141)

Now a cool “quiver” plot

if flag:
    matplotlib.pyplot.rcParams["figure.figsize"] = [10.00, 8.00]
    matplotlib.pyplot.rcParams["figure.autolayout"] = True
    im = matplotlib.pyplot.imread("Fig5.18-Coordinates.png") # base image
fig, ax = matplotlib.pyplot.subplots()
if flag: 
    im = ax.imshow(im, extent=[-25, 1000, -125, 800])# sets X and Y plot window of basemap
#fig.set_size_inches(14, 7)
levels=[5.50,5.75,6.00,6.25,6.50,6.75]
CS = ax.quiver(XX[::6, ::6], YY[::6, ::6], UU[::6, ::6], VV[::6, ::6], units='width')
#ax.clabel(CS, inline=2, fontsize=12)
ax.set_title('Vector Plot from Gridded Data File')
ax.set_xlim([0,1000])
ax.set_ylim([0,800]);
ax.plot([xxl,xxh],[yyl,yyl],color="black")
ax.plot([xxl,xxh],[yyh,yyh],color="black")
ax.plot([xxl,xxl],[yyl,yyh],color="black")
ax.plot([xxh,xxh],[yyl,yyh],color="black");

#matplotlib.pyplot.quiver(XX[::6, ::6], YY[::6, ::6], UU[::6, ::6], VV[::6, ::6], units='width')
#matplotlib.pyplot.streamplot(XX, YY, UU, VV, density=0.5, linewidth=2, color=None)
../../_images/8367502528ff7611e622a7cd79e76576fc5b3ff90bb1824f2b51e8dba9656e79.png 
if flag:
    matplotlib.pyplot.rcParams["figure.figsize"] = [10.00, 8.00]
    matplotlib.pyplot.rcParams["figure.autolayout"] = True
    im = matplotlib.pyplot.imread("Fig5.18-Coordinates.png") # base image
fig, ax = matplotlib.pyplot.subplots()
if flag: 
    im = ax.imshow(im, extent=[0, 1000, -125, 800])# sets X and Y plot window of basemap
#fig.set_size_inches(14, 7)
levels=[5.50,5.75,6.00,6.25,6.50,6.75]
CS = ax.streamplot(XX, YY, UU, VV,  density=0.5, linewidth=1, color=None, arrowsize=3)
#ax.clabel(CS, inline=2, fontsize=12)
ax.set_title('Vector Plot from Gridded Data File')
ax.set_xlim([0,1000])
ax.set_ylim([0,800]);
ax.plot([xxl,xxh],[yyl,yyl],color="black")
ax.plot([xxl,xxh],[yyh,yyh],color="black")
ax.plot([xxl,xxl],[yyl,yyh],color="black")
ax.plot([xxh,xxh],[yyl,yyh],color="black");
../../_images/6f20600fd5392667e2d947a2bc5a87216d5bd74266b7d0ee47af9901b12d541c.png

Now we can implement a particle tracking script. Using the spreadsheet as a guide we can write it in python, and access the cool graphics.

Here is the basic script, one simply has to wrap it into a time-stepping loop to produce a trajectory. Multiple particles can be managed

# U array x-velocity at XG,YG
# V array y-velocity at XG,YG
# XG array X-value of cell center
# YG array Y-value of cell center
# XP X-value of particle position 
# YP Y-value of particle position
# UP x-velocity of particle 
# VP y-velocity of particle
# TX x-component particle trajectory
# TY y-component particle trajectory

import math
verbose=False
terse=False
deltaT = 100
etime = 0
numTime = 42
XP = []
YP = []
UP = []
VP = []
TX = [] #trajectory vector
TY = []
np = 1 # Total particles 
XP.append(900)
YP.append(550)
UP.append(0)
VP.append(0)
ip=np-1 
print("                   Initial Particle Position",round(XP[ip],2),round(YP[ip],2),round(UP[ip],2),round(VP[ip],2),round(etime,2))

# move particles this time step
for it in range(numTime):
    for ip in range(np):
    # Build Particle Distance Table
        dist = []
        index = []
        count = 0
        for j in range(ncols):
            for i in range(nrows): #range(nrows)
                dist.append(math.sqrt((XX[i][j]-XP[ip])**2 + (YY[i][j]-YP[ip])**2))
                index.append([count,i,j]) # use to find i,j for a given index
                count = count +1
    # find closest cell
        for i in range(count):
            if dist[i] <= min(dist):
                #print(index[i],dist[i])
                ixx=index[i][1]
                jyy=index[i][2]
    # use nearest cell assignment - aka simple scheme
                UP[ip] = UU[ixx][jyy]
                VP[ip] = VV[ixx][jyy]
                if verbose: print("Particle Position and Velocities Before Move",round(XP[ip],2),round(YP[ip],2),round(UP[ip],2),round(VP[ip],2),round(etime,2))
                break #exits the loop - we stop at the first nearest cell encountered
    # move the particle
        XP[ip]=XP[ip]+UP[ip]*deltaT
        YP[ip]=YP[ip]+VP[ip]*deltaT
        etime=etime+deltaT
        if terse: print(" Particle Position and Velocities After Move",round(XP[ip],2),round(YP[ip],2),round(UP[ip],2),round(VP[ip],2),round(etime,2))
        TX.append([ip,XP[ip],etime])
        TY.append([ip,YP[ip],etime])

                   Initial Particle Position 900 550 0 0 0

if flag:
    matplotlib.pyplot.rcParams["figure.figsize"] = [10.00, 8.00]
    matplotlib.pyplot.rcParams["figure.autolayout"] = True
    im = matplotlib.pyplot.imread("Fig5.18-Coordinates.png") # base image
fig, ax = matplotlib.pyplot.subplots()
if flag: 
    im = ax.imshow(im, extent=[0, 1000, -125, 800])# sets X and Y plot window of basemap
#fig.set_size_inches(14, 7)
levels=[5.50,5.75,6.00,6.25,6.50,6.75]
CS = ax.streamplot(XX, YY, UU, VV,  density=0.5, linewidth=1, color=None, arrowsize=3)
#ax.clabel(CS, inline=2, fontsize=12)
ax.set_title('Vector Plot from Gridded Data File\n' +\
            'Each Marker is ' + str(deltaT) + ' days\n' +\
            'Total Time is ' + str(etime) + ' days')
ax.set_xlim([0,1000])
ax.set_ylim([0,800]);
ax.plot([xxl,xxh],[yyl,yyl],color="black")
ax.plot([xxl,xxh],[yyh,yyh],color="black")
ax.plot([xxl,xxl],[yyl,yyh],color="black")
ax.plot([xxh,xxh],[yyl,yyh],color="black")
xtrajectory = [sublist[1] for sublist in TX]
ytrajectory = [sublist[1] for sublist in TY]
ax.plot(xtrajectory[0],ytrajectory[0],marker="o",color="blue",markersize=24)
ax.plot(xtrajectory,ytrajectory,marker="o",color="red")
;

''
../../_images/a0771ff432ba6d92891fa5fbd8eb429aa4cb535c71fe44c47cb317095512c1a9.png

Extensions#

  1. One could certainly add a lot more particles.

  2. If the parcels represent some quantify of mass, concentrations can be directly approximated.

The method is not limited to groundwater hydrology, and an example of another use is in the references.

References#

  1. Cleveland, T.G., 1991. A Comparison of Sampling Design Criterion using a Lagrangian Particle Tracking Model for Transport in Porous Media, Final Report to Houston Advanced Research Center, Woodlands, Texas. 54p.

  2. Ahlstrom, S W, Foote, H P, Arnett, R C, Cole, C R, and Serne, R J. Multicomponent mass transport model: theory and numerical implementation (discrete-parcel-random-walk version). United States: N. p., 1977. Web. doi:10.2172/7083383.

  3. Modeling plume behavior for nonlinearly sorbing solutes in saturated homogeneous porous media in Advances in Water Resources 21 (1998) 487-498 Elsevier Science Limited

  4. Bear, J. (1972) Dynamics of Fluids in Porous Media McGraw Hill (pp. 236-241)

  5. Cleveland, T.G., Thompson, D.B., Fang,X., He, X. 2007. Synthesis of Unit Hydrographs from a Digital Elevation Model. Urban Runoff Modeling: Intelligent Modeling to Improve Stormwater Management, ECI Specialty Conference, July 22-27, Humboldt State University, Arcata, California.

Contents