Lecture 10: Finite Difference Methods

Numerical solution to the Laplace equation

The above simple situation almost never exists in reality. When real situations are incoroporated in the equation, the equation becomes impossible to solve. We then can only approximate solutions.

Numerical methods are mathematical techniques to approximate solutions.

Finite difference is a commonly applied numerical method.

The basic concept is that we break the continuous domain into a series of finite domains that approximate the continuous curve with a line (mathematicians would cringe at this explanation). Consider a power function…

Consider our 1D groundwater flow problem. Whereas our analytical solution allows us to solve for h at any x, the finite difference method requires us to discretize the problem.

Derive the finite difference approximation for the 1D Laplace equation…

hi = (hi-1+hi+1)/2

Discuss constraints and assumptions

-linear potentiometric surface

-What does this mean for grid size?

Demonstrate 1D Excel solution

-Introduce

-Iteration

-Guessing and adjusting to solve simultaneous equations

-Convergence

Numerical solution to the Laplace equation using Microsoft Excel.

The purpose of steady-state groundwater modeling is to simulate the distribution of head in a flow system from just a few points where head is known. One of the simplest systems, and the one we will focus on, is a 2-d confined, homogeneous, isotropic aquifer. For this situation, the governing equation is

Recall that the analytical solution to the 1-D version of this equation is simply a linear equation

h = ax + b

where a and b are constants that are evaluated by plugging in boundary conditions (locations where h is known).

In 2-D, the analytical solution is more complicated and analytical solutions exist only for the simplest cases. Consider the Toth problem – a flow system with 3 impermeable boundaries with a linearly sloping potentiometric surface:

The analytical solution to the Toth problem is shown on page 19 of Wang and Anderson. You notice that there is a large jump in complexity from 1d to 2d.

Instead, we replace the continuous domain with a finite grid.

See page 22 of Wang and Anderson…

Hi,j = (1/4)(hi-1,j + hi+1,j + hi,j-1 + hi,j-1)

The subscript i and j refer to the x and y directions, respectively. The equation simply says that the value of cell i,j is simply the average of the four surrounding cells.

The above equation tells us the relationship of head in a cell to the other cells, but how do we get the absolute values?

…boundary conditions. We have to know the actual value of head somewhere in the problem. Then the above equation will ripple through the problem until it converges.

Bounday conditions can be;

Known head – called Dirichlet condition
Known flux – called Neumann condition

The water table is an example of a known head boundary condition. A no-flow boundary is an example of a known flux boundary.

Let’s work through an example.

Example 1.

Consider a cross-section of a 2-d groundwater flow field 1000 ft long and 1000 feet deep. The aquifer is unconfined, homogeneous, and isotropic. The elevation of the water table is 1000 feet at the west end of the problem and 1100 at the east end of the problem. The two vertical boundaries and the bottom boundary are no-flow barriers. Assume that the water table has a linear slope and that the head in the aquifer is desribed by the LaPlace equation. Use a finite-difference approximation of the Laplace equation to construct a map of the flow field with equipotential lines and flowlines. Hand in your map electronically.

1. Start Microsoft Excel

2. Set up the grid

a. Enter cell dimensions along the top and left sides of the problem domain.

b. Adjust cell widths and heights so that cells are square (isotropic).

c. Highlight the cell borders. These cells are outside the problem domain and are here just for scale reference.

d. Format interior cells so that the contents appear in the center of the cell.

3. Shade the boundary cells a color of your choice for the time being.

4. Enter finite difference equations in all interior cells. Note that 0’s appear in all cells. That’s because there are no boundary conditions to provide a seed value. Just for illustration enter any number in the top cell of the second column. Press F9 to activate the equations. You will see the solution ripple through all cells once. The problem is not done, however, because when a cell value changes it changes the values of the surrounding cells which changes itself… Hold down the F9 key continuously to see what I mean. The solution will slow down and eventually stop when the convergence criteria is met.

5. Set the convergence criteria. Click Tools, Options, Calculation. Check the iteration box if it is not check. This will allow for circular references. The equations will be resolved until the change in values between iterations in every cell is less than a specified value. Enter a value of 0.001 in the maximum change box.

6. Get rid of the bogus value you enter in step 4. The interior cells will still have numbers. We will use these as initial conditions.

7. Enter boundary conditions.

a. Enter 1000 in the upper left corner and 1100 in the upper right corner.

b. Fill in the equation for a line for all other upper cells.

Y = mx+b

c. For no-flow boundaries around perimeter (excluding corners)

hx,0 = (hx-1,y + 2*hx,y+1 + hx+1,y)/4

h0,y = (hx,y+1 + 2*hx+1,y + hx,y-1)/4

hx,y = (hx,y-1 + 2*hx-1,y + hx,y+1)/4

d. For bottom corners

hx,y = (2*hy+1,x + 2*hx+1,y)/4 or hx,y = (2*hy+1,x + 2*hx-1,y)/4

8. Hold down F9 key until there are no more changes in any cell.

9. Contour the results.

10. Draw flowlines

Problem 1.

Consider a plan view of a confined, isotropic, homogeneous aquifer that is 1000 m wide and 1000 m long. The west boundary has a constant head of 1500 m and the east boundary has a constant head of 1700 m.

a. Construct a potentiometric map and flow lines for the condition described. If the hydraulic conductivity of the aquifer is 10-2cm/s and the thickness is 200 m, calculate the flow rate through the aquifer. If the porosity is 0.3, calculate the average linear velocity of groundwater.

b. Assume a well is installed in the center cell and after equilibrium it maintains a constant head of 1450. Construct the new potentiometric map and flow lines. Hand in your map electronically.

Problem 2. Consider a cross-section of a 2-d groundwater flow field 1000 ft long and 1000 feet deep. The aquifer is unconfined, homogeneous, and isotropic. The elevation of the water table is 1000 feet at the both ends and 1100 in the center. Assume the water table has a linear slope from the center to the ends. The two vertical boundaries and the bottom boundary are no-flow barriers. Assume that the water table has a linear slope and that the head in the aquifer is desribed by the LaPlace equation. Use a finite-difference approximation of the Laplace equation to construct a map of the flow field.

Problem 3.

Repeat example 1, but simulate the influence of undulating topography on the regional sloping water table. Assume that the water table mimics the topography and can be approximated with a sin wave governed by

hx = (z0 + Bx/L + bsin(2*pi*x/n))

where B is the total relief of the water table

L is the total length in the x direction of the domain

z0 is the elevation of the water table at a known point

x is the distance from z0

b is the amplitude of the sin wave

n is the number of oscillations in the sin wave (L divided by the number of flow cells)