Solving 2-D Groundwater Head: In Unsteady State Condition Of Confined Aquifer Using C#

Introduction

This code is useful in the field of Environmental and Hydrogeology Sciences. Where most groundwater problems appear in most cases such as calculating the groundwater budget under the stresses of discharge/recharge (pumped or injected wells).

Background

Finite difference approximation is a fundamental numerical math method that helps to solve the problem of complex equations that function natural physical problem. Especially in the problem of transient condition (time dependent function), and there are three general widely used finite difference schemes.

1. Explicit Differentiation

h(i,j,k+1) = rh(i-1,j,k)+(1-4r)h(i,j,k)+rh(i+1,j,k)+rh(i,j-1,k)+rh(i,j+1,k)

2. Implicit

h(i,j,k) = rh(i-1,j,k+1)+(1+4r)h(i,j,k+1)+rh(i+1,j,k+1)+rh(i,j-1,k+1)+rh(i,j+1,k+1)

3. Crank Nicholson

rh(i-1,j,k)+(2-4r)h(i,j,k)+rh(i+1,j,k)+rh(i,j-1,k)+rh(i,j+1,k) = 
               rh(i-1,j,k+1)+(2+4r)h(i,j,k+1)+rh(i+1,j,k+1)+rh(i,j-1,k+1)+rh(i,j+1,k+1)

And finally, combining all three methods by adding Lamda coefficient (code below):

h(i,j,k+1) = [h(i,j,k)+r*lamda*[h(i-1,j,k)-4h(i,j,k)+h(i+1,j,k)+
h(i,j-1,k)+h(i,j+1,k)]+r(1-lamda)*[h(i-1,j,k+1)+h(i+1,j,k+1)+h(i,j-1,k+1)+h(i,j+1,k+1)]]/(1+4r(1-lamda)

where:

• r = D*DT/DX^2
• D = T/S
• T: is the transmissivity of the aquifer (m2/day)
• S: aquifer sorage coefficient
• DX: change in distance in X-direction
• DT: time increment
• lamda:
  • if = 0 ----> Implicit, if = 0.5 ----> Crank Nicholson, if = 1 ----> Explicit

Using the Code

In this code, we will solve the groundwater head of confined aquifer distribution in which a sudden drop occurs to groundwater at initial condition was 8.8 m distributed in all the aquifer domain. The water flow from upstream boundary with constant head of 8.8 m towards the downstream boundary 3.2 m.

The figure below illustrates the conceptual model of the aquifer:

In the code below, Gauss Seidel Iteration method was applied to solve the finite difference scheme.

            double T = 2;            // Transmissivity (m2/day)
            double S = 0.0003;       // Storativity 
            double DT = 0.5;         // Time Step Increment (day)
            double DX = 100;         // Spatial Increment (m)
            double lamda = 0;        // if lamda = 0 ---> Implicit , 
                                     // if lamda = 0.5 ---> Crank Nicolson, if lamda = 1 ---> Explicit 
            double D = T / S;
            double H0 = 8.8;         // Initial Condition (m)
            double H1 = 8.8;         // Upstream Head Boundary (m)
            double H2 = 3.2;         // Downstream Head Boundary (m)
            int IT = 100;            // Number of Iterations
            double ERR = 1;          // Error Function
            double TOL = 0.001;      // Tolerance 
            double r = (D * DT) / System.Math.Pow(DX, 2);

            double[, ,] H = new double[10, 10, 18];

            for (int i = 0; i <= 9; i++)
            {
                for (int j = 0; j <= 9; j++)
                {
                    int k = 0;

                    H[i, j, k] = H0;
                }
            }

            for (int i = 0; i <= 9; i++)
            {
                int j = 0;

                for (int k = 1; k <= 17; k++)
                {

                    H[i, j, k] = H1;
                }
            }

            for (int i = 0; i <= 9; i++)
            {
                int j = 9;

                for (int k = 1; k <= 17; k++)
                {

                    H[i, j, k] = H2;
                }
            }

            while (ERR > TOL)

                for (int n = 0; n < IT; ++n)
                {
                    for (int i = 1; i <= 8; i++)
                    {
                        for (int j = 1; j <= 8; j++)
                        {
                            for (int k = 0; k <= 16; k++)
                            {
                                double OLD = H[i, j, k + 1];
                                H[i, j, k + 1] = (H[i, j, k] + (r * lamda) * (H[i - 1, j, k] - 
                                (4 * H[i, j, k]) + 
                                H[i + 1, j, k] + H[i, j - 1, k] + H[i, j + 1, k]) +
                                (r * (1 - lamda)) * (H[i - 1, j, k + 1] + H[i + 1, j, k + 1] + 
                                H[i, j - 1, k + 1] + 
                                H[i, j + 1, k + 1])) / (1 + (4 * r * (1 - lamda)));
                                H[0, j, k + 1] = H[2, j, k + 1];
                                H[9, j, k + 1] = H[7, j, k + 1];
                                ERR = System.Math.Abs(H[i, j, k + 1] - OLD);
                            }
                        }
                    }
                }

            Console.WriteLine("Groundwater Head values in (m)");
            Console.WriteLine(IT.ToString());
            Console.WriteLine(ERR.ToString());
            for (int i = 0; i <= 9; i++)
            {
                for (int j = 0; j <= 9; j++)
                {
                    int k = 17;

                    Console.Write(string.Format("{0} ", H[i, j, k]));
                }
                Console.Write(Environment.NewLine + Environment.NewLine);
            }
            Console.ReadLine();
        }
    }
}

Output

Head distribution at the end of time (t = 18), with number of iterations = 100

tolerance = 0.002

Points of Interest

This code was implemented with C# language using console application to solve the PDE in the environmental sciences with the application of FDM in two-dimensions.

References

• Mary P Anderson. Introduction to Groundwater Modeling: Finite Difference Finite Element
• Amr El-Feki. Prof in Environmental Sciences K A U

History

• 4^th October, 2018: Initial version