8. Modeling Advection#
Readings#
Bear, J. (1972) Dynamics of Fluids in Porous Media McGraw Hill (pp. 236-241)
Koutitas,C.G. (1983) Elements of computational hydraulics Pentech Press ; New York : Distributed in the USA by Chapman and Hall Link is to internet book archive.
Modeling advective transport in groundwater is a critical component of hydrogeology and environmental science, as it allows us to understand and predict the movement of contaminants and the transport of solutes through subsurface aquifers. Advective transport refers to the bulk movement of water and the substances it carries, driven primarily by hydraulic gradients and pressure differences within the groundwater system. This process plays a fundamental role in determining the fate and transport of pollutants, nutrients, and other substances in the subsurface environment.
A variety of widely used methods and techniques are used; These methodologies are essential for comprehending the dynamics of groundwater flow and the transport of substances within aquifers. Two popular methods for modeling advective transport in groundwater are numerical modeling and analytical solutions.
Numerical modeling stands as a cornerstone in the field, with models like MODFLOW (MODular Finite-Difference Ground-Water Flow Model) at the forefront. Numerical models divide the subsurface into discrete grid cells, considering various factors such as hydraulic conductivity, boundary conditions, and aquifer heterogeneity. These models utilize numerical algorithms to simulate the behavior of groundwater flow and contaminant transport over time. By incorporating data from field investigations, they enable scientists to make accurate predictions regarding the movement and fate of contaminants in complex geological settings. Numerical models are exceptionally versatile and can handle a wide range of scenarios, making them an indispensable tool for hydrogeologists and environmental scientists.
On the other hand, analytical solutions offer a more simplified approach for specific scenarios. These solutions are often employed when dealing with idealized geometries and boundary conditions. For instance, the advection-dispersion equation, a fundamental equation in groundwater transport, can be solved analytically for simple scenarios. Analytical solutions provide quick insights into transport phenomena and are valuable for initial assessments and educational purposes. However, they may not capture the complexities of real-world aquifers as effectively as numerical models.
Methods#
A reductionist point of view lets us categorize practical modeling methods into two categories: particle (parcel) tracking and eulerian finite- difference,element, or volume methods (marker-in-cell is a hybrid). The Eulerian methods almost always employ some kind of flux limiter to maintain proper mass balances, upwinding is one such approach. Particle tracking and upwind finite difference approaches are two distinct methods for modeling advective transport in groundwater, each with its own set of advantages and limitations. Let’s compare and contrast these approaches:
Particle Tracking:#
Conceptual Approach: Particle tracking models simulate the movement of individual particles within the groundwater flow field. This approach provides a more intuitive understanding of how contaminants or solutes move through the subsurface.
Accuracy in Complex Flow: Particle tracking is well-suited for modeling complex groundwater flow scenarios, such as in fractured aquifers or regions with varying flow directions. It can accurately represent irregular flow patterns.
Limitations: However, particle tracking may become computationally intensive when modeling a large number of particles over extended periods, making it less practical for long-term simulations.
Adaptive Resolution: It allows for adaptive grid refinement, where particles can be concentrated in regions of interest, offering higher resolution where needed.
Lagrangian Approach: It follows the Lagrangian approach, where particles move with the flow, which can be advantageous for understanding transport phenomena from a particle’s perspective.
Upwind Finite Difference:#
Grid-Based Approach: Upwind finite difference methods are grid-based, where the domain is divided into grid cells, and equations are solved at discrete locations within each cell.
Computational Efficiency: Upwind finite difference methods are computationally efficient for large-scale models and long-term simulations. They are widely used for regional-scale groundwater flow and transport modeling.
Stability: They are known for their numerical stability, which ensures that the solution remains accurate even when using relatively large time steps.
Accuracy in Uniform Flow: Upwind finite difference methods may perform very well in uniform flow conditions, but they can encounter challenges in accurately representing complex flow patterns, especially in situations where there are rapid changes in flow directions.
Eulerian Approach: These methods follow an Eulerian approach, where flow and transport equations are solved at fixed grid locations, making them less intuitive from a particle’s perspective.
In summary, the choice between particle tracking and upwind finite difference approaches for modeling advective transport in groundwater depends on the specific characteristics of the groundwater system and the research objectives:
Particle tracking is preferred for understanding transport at a detailed, particle-level scale and for modeling complex flow patterns.
Upwind finite difference methods are suitable for large-scale, long-term simulations where computational efficiency and stability are essential, especially in cases of relatively uniform flow.
Often, a combination of both methods may be employed in a modeling study to take advantage of their respective strengths, using particle tracking for detailed insights in specific areas and upwind finite difference methods for broader-scale simulations.
The remainder of this notebook presents simple examples of each method.
Note
While the methods in this notebook would be fine, practical “for money” modeling will usually employ specialized software; especially if lawyers will somehow get involved.
To build models, we simply apply the Reynolds Transport Theorem to the advective flux component to produce something like:
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
INDEXfunction allows us to select the correct values of velocity if we know which cell the particle resides in.The
MATCHfunction allows us to compare values in an array and determine the position in the array that these values are found.
Thus theMATCHfunction 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#
Cleveland, T.G. (1997) Particle Tracking Spreadsheet (10X10) (Excel)
Cleveland, T.G. (1997) Particle Tracking Spreadsheet (20X20) (Excel)
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:
Which well is expected to be the most contaminated.
The groundwater velocity and seepage velocity across the plume.
The duration that the source has been contaminating the aquifer (neglect dispersion, diffusiom, and adsorption).
The flow rate across the plume.
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.
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
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)
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('Streamline 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");
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(290)
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 290 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('Streamline 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")
;
''
Extensions#
One could certainly add a lot more particles.
If the parcels represent some quantify of mass, concentrations can be directly approximated.
The method is not limited to groundwater hydrology, it is useful for surface water hydrology to find flow paths for time-of-concentration calculations; extension to the 3rd spatial dimension is also straightforward.
Eulerian (Finite-Difference) Approach#
The transport equations are partial differential equations (PDEs), which describe the evolution of a quantity over space and time. One of the most common methods for solving these equations numerically is the finite-difference method (FDM). This approach discretizes space and time into a grid and approximates the derivatives of the governing equations.
Finite-Difference Method#
In the finite-difference method, spatial and temporal domains are discretized into small steps. The derivatives in the governing transport PDE are approximated using the differences between neighboring grid points.
For example, consider the 1D advection equation (a typical transport equation): $\(\frac{\partial C}{\partial t} = - V \frac{\partial C}{\partial x}\)$
where \(C(x,t)\) is the transported quantity (e.g., concentration), \(V\) is the advection velocity, and \(t\) and \(x\) are time and space, respectively.
To solve approximate the solution using finite-difference, we could approximate the derivatives as:
However, the simple central difference approximation for the spatial derivative may lead to non-physical oscillations or instability, particularly when solving hyperbolic dominanted transport equations like the advection equation.
Upwinding in Hyperbolic-Dominant PDEs#
In transport modeling, hyperbolic PDEs (like the advection equation) describe wave propagation and are characterized by the dominance of the advection term. When solving such equations numerically, the choice of difference scheme is critical for stability and accuracy.
Upwind differencing is a technique specifically designed to handle hyperbolic PDEs. The idea is to use information from the direction of the flow (upstream) to approximate the spatial derivative. This prevents the introduction of non-physical oscillations by ensuring that the finite-difference scheme respects the direction of the characteristic waves (i.e., the direction in which information propagates). Upwind Scheme
For a simple 1D advection equation, if the flow velocity \(V>0\) (information travels to the right), we use the backward difference to approximate the spatial derivative: $\(\frac{\partial C}{\partial x} \approx \frac{C_{i}^{n} - C_{i-1}^n}{\Delta x}~\text{(spatial derivative, backward difference)}\)$
If \(V<0\) (information travels to the left), we use the forward difference:
This technique, called upwinding, introduces numerical dissipation that stabilizes the solution and prevents oscillations, but it does so at the cost of some accuracy (it can be less accurate for smooth solutions compared to central differencing).
Why Upwinding is Important#
Upwinding is essential in transport modeling because most real-world problems involve advection-dominated transport, where sharp fronts or discontinuities are common (e.g., contaminant plumes or shock waves). Without upwinding, traditional finite-difference methods can result in:
Non-physical oscillations: especially when trying to capture sharp gradients.
Instability: numerical solutions can blow up if the scheme is not carefully designed for hyperbolic systems.
Artificial diffusion: even stable solutions will have apparent added diffusion that is absent from the physical system.
By incorporating upwinding, we align the numerical method with the physics of the problem, ensuring that information propagates in the correct direction, thereby stabilizing the simulation.
A simple example follows to illustrate the effect of upwinding and no upwinding.
Upwind formulation#
A simple upwind formulation for 1D system is
The parts with the absolute value are a sneaky way to handle the upwinding automatically.
As long as \(V_i <= \frac{\Delta t}{\Delta x}\) this difference scheme is stable, it is also reasonably accurate if \(V_i \approx \frac{\Delta t}{\Delta x}\)
Below is a simple spreadsheet implementation of the explicit scheme above.
In MODFLOW6 transport model, a choice of upwind (called upstream in the toolkit) tries to match internal time steps to maintain numerical stability. As you will see in class when we run a spreadsheet example. Its impossible for perfect accuracy and stability, but it is possible to stay close. We will see examples of numerical diffusion as we depart from the accuracy requirement.
Centered formulation#
This scheme will exhibit substantial numerical diffusion in part because it is a non-upwind formulation.
Below is a simple spreadsheet implementation of the explicit, centered in space scheme above.
The spreadsheet above nicely exhibits numerical diffusion, apparent in the profile plot as well as oscillation in the time history. Professional tools address this phenomenon by a variety of approaches, upwinding and flux limiters being common reliable methods, a method of characteristics is another way to deal with the issue.
Note
Need to build python examples
Numerical Diffusion When solving advection-dominated partial differential equations (PDEs) using numerical methods, one common issue is numerical diffusion (also called numerical dissipation). This artifact of the numerical scheme can significantly affect the accuracy of the solution, particularly when sharp gradients or fronts are present, as in contaminant transport or shock waves.
Numerical diffusion is especially problematic in non-upwind formulations, such as those using central differencing for spatial derivatives. To understand why this happens, let’s review what numerical diffusion is and why it arises in these schemes. What is Numerical Diffusion?
Numerical diffusion refers to the artificial smoothing or blurring of a solution in a numerical model, which leads to a loss of sharpness in the advective transport of quantities like contaminants, pollutants, or heat. This effect is not physical but arises from the discretization of the advection equation.
The equation should transport \(C(x,t)\) without changing its waveform (if no other processes, such as diffusion or dispersion, are acting). However, many numerical schemes, particularly non-upwind methods, introduce spurious diffusion into the solution.
Central Differencing and Numerical Diffusion
One popular method for approximating the spatial derivative is central differencing which is used in the scheme just preceding. While central differencing has second-order accuracy, meaning it can provide higher accuracy for smooth solutions, it is unstable when applied to advection-dominated problems without additional stabilization mechanisms. Specifically, central differencing often leads to oscillations and numerical diffusion in the solution when dealing with sharp fronts or discontinuities.
Why Numerical Diffusion Occurs
Numerical diffusion arises from the fact that non-upwind methods, like central differencing, do not adequately account for the direction of the flow in advection-dominated problems. In these cases, the finite-difference scheme inadvertently mimics a diffusive process. This can be understood as follows:
Information propagation: In hyperbolic PDEs like the advection equation, information propagates in a specific direction (determined by the sign of cc). Non-upwind methods like central differencing treat the spatial derivative symmetrically, which leads to a lack of proper alignment with the flow direction.
Smoothing effect: When central differencing is used in such cases, it smooths the solution because it averages information from both directions (upstream and downstream) of the grid point. This acts similarly to a diffusion process, where sharp gradients are smeared over time.
Impact of Numerical Diffusion
Numerical diffusion is particularly problematic for sharp fronts, such as contaminant plumes, where it can cause:
Excessive smoothing of the front: The sharp boundary of the plume becomes diffused over a larger area than physically expected.
Loss of mass conservation: In some cases, numerical diffusion can lead to non-conservative solutions, where the total amount of the transported quantity is not accurately preserved.
For example, if you are modeling the transport of a pollutant in groundwater using central differencing, the sharp concentration gradient between the polluted and unpolluted regions can be artificially blurred. This causes the plume to spread out more than it physically should.
Mitigating Numerical Diffusion
To reduce or eliminate numerical diffusion, upwinding schemes are commonly used. The upwind finite-difference method ensures that the spatial derivative is calculated using values from the upstream direction, thus preventing the unwanted smoothing that central differencing introduces. Although upwinding itself introduces some numerical dissipation, this dissipation is somewhat controlled and aligned with the flow direction, making it more suitable for advection-dominated problems.
Marker-in-Cell and Method of Characteristic models are hybrids that leverage a particle tracking approach to locate mean values, then dispersion about that mean value to honor the hyperbolic component while preserving the diffusion-like nature of the non-advective terms.