Suppose that the velocity field is given (e.g. it may have been estimated from water level measurements or simulated with a groundwater flow model).  We wish to write a MATLAB program to keep follow a number of particles as they move from left to right.

In this problem we need to consider variability over space as well as time.  This requires the definition of a rectangular spatial grid such as the one shown here.  Velocities are defined at the intersections of the grid lines (or nodes).  The well is located at the origin (x = 0, y = 0).  Particles are released at specified (x, y) coordinates at the initial time (t = 0).  As the particles move through the velocity field their x and y coordinates change.  These coordinates need to be recomputed at each time step.  The particles can be visualized by plotting  their coordinates as discrete points on an (x,y) graph which spans the computational grid.  Particles that leave the grid, either across the sides or through the well, must be removed.  Our objective is to plot the particles in a series of snapshots taken at different times so we can see how they move towards or past the well. 

. 

Two types of calculations are required for this problem: 1) calculation of the velocity at each grid node and 2) calculation of the x and y coordinates of each particle at each time.  The velocity at each node is the sum of two components: 1) a regional velocity u, which we assume is a constant pointing in the x direction and 2) the well pumping velocity, which points radially towards the well.  The well pumping velocity component may be derived by applying the principle of mass conservation.  Suppose that the grid is viewed as an array, with i the row (or grid line aligned with the x direction) and j the column (or grid line aligned in the y direction).  The resulting equations for the x and y  velocities at node (i,j) are:

Once the x and y velocities are computed, the coordinates of each particle at the new time (t+1) are related to the coordinates at the old time (t) as follows:

Translate the solution into MATLAB code:  The downloadable MATLAB code tracking.m solves the stated problem:

 

function tracking % Program to simulate particle movement in the capture zone of a % pumping well % Define velocity grid nx=15 ny=15 u=1.5  % regional groundwater velocity (constant in x direction) Q=-20. % well pumping rate dx=1.  % x grid spacing dy=1.  % y grid spacing d=1.   % flow depth % define grid-related arrays x=dx*[-(nx-1)/2:(nx-1)/2]; y=dy*[(ny-1)/2:-1:-(ny-1)/2]; xmin=x(1); xmax=x(nx); ymax=y(1); ymin=y(ny); % bigx and bigy define x and y corrdinates respectively, of every % node in grid bigx=ones(ny,1)*x; bigy=y'*ones(1,nx); % compute velocities (from mass balance) denom=2*pi*d*sqrt(bigx.^2+bigy.^2); % radial velocity denom. zeroflag=(denom==0); % check for divide by zero denom(zeroflag)=1.; % set zero elements of denom =1. vr=Q./denom;  % compute vr  vr(zeroflag)=0.; % reset vr = 0 where denom was originally zero vx=u*ones(ny,nx)+vr.*cos(atan2(bigy,bigx)); vy=vr.*sin(atan2(bigy,bigx)); % plot velocity field close all figure quiver(bigx,bigy,vx,vy) axis([xmin, xmax, ymin, ymax]) % specify number of paricles and initial locations % y locations distributed randomly np=100 yminp=ymin+0.10*(ymax-ymin); ymaxp=ymin+0.90*(ymax-ymin); xp(1:np)=x(2); yp(1:np)=yminp*ones(1,np)+(ymaxp-yminp)*rand(1,np); % specify number and size of tracking time steps nt=40 dt=0.5 % start time loop for t=1:nt    k=1;    while (k<=np)    % interpolate velocities from grid to particle locations       vxint=interp2(bigx,bigy,vx,xp(k),yp(k));       vyint=interp2(bigx,bigy,vy,xp(k),yp(k));       % compute nominal new locations       ytemp=yp(k)+vyint*dt;       xtemp=xp(k)+vxint*dt;       % move particles while they remain within grid       if(((xtemp~=0.)|(ytemp~=0.))&(xtemp>=xmin)&(xtemp<=xmax)...             &(ytemp>=ymin)&(ytemp<=ymax))          xp(k)=xtemp;          yp(k)=ytemp;       % remove particles that have left domain       else          xp(k)=[];          yp(k)=[];          np=np-1       end       k=k+1;    end    % plot particle locations at selected times    if((t==1)|(mod(t,2)==0))       figure       plot(xp,yp,'o')       axis([xmin, xmax, ymin, ymax])    else    end end return