Fundamental Approach in Groundwater Flow and Solute Transport Modelling Using the Finite Difference Method

Water below the ground surface is usually referred to as groundwater. Groundwater has recently become a major source of water supply in almost every sector. Over-dependence on it for many purposes has led to its’ over-exploitation, and this has led to much concern for groundwater assessment and management. For a proper assessment and management of groundwater resources, a thorough understanding of the complexity of its processes is quite essential. Expansion of human activities causes dispersion of pollutants in the subsurface environment. The fate and movement of dissolved substances in soils and groundwater has generated considerable interest out of concern for the quality of the subsurface environment Groundwater flow and transport analysis have been an important research topic in the last three decades. This is as a result of many geo-environmental engineering problems having direct or indirect impact by groundwater flow and solute transport. Solute transport by flowing water (dissolved suspended particles) has broad impact in environmental protection and resource utilization via groundwater contamination. The leaching (displacement) of salts and nutrients in soils also has an impact on agricultural production. To predict the contaminant migration in the geological formation more accurately, many analytical solutions for partial differential equations exist but because of the difficulty in obtaining analytical solutions, numerical solutions are more generally used. Numerical solutions are often more difficult to verify, so mathematical model error has to be kept as small as possible. Many techniques for solving numerically the solute transport (advection-dispersion) equation are used such as the Finite Element Method (FEM), Finite Difference Method (FDM), Boundary Element Method (BEM), Fuzzy Sets Approach (FSA), Artificial Neural Networks method (ANN), Particle Tracking Method (PTM), Random Walk Method(RWM), Integrated Finite Difference Method (IFDM) etc. Conventional analysis of groundwater flow is generally made by using the relevant physical principles of Darcy's law and mass balance.


Introduction
Water below the ground surface is usually referred to as groundwater. Groundwater has recently become a major source of water supply in almost every sector. Over-dependence on it for many purposes has led to its' over-exploitation, and this has led to much concern for groundwater assessment and management. For a proper assessment and management of groundwater resources, a thorough understanding of the complexity of its processes is quite essential. Expansion of human activities causes dispersion of pollutants in the subsurface environment. The fate and movement of dissolved substances in soils and groundwater has generated considerable interest out of concern for the quality of the subsurface environment Groundwater flow and transport analysis have been an important research topic in the last three decades. This is as a result of many geo-environmental engineering problems having direct or indirect impact by groundwater flow and solute transport. Solute transport by flowing water (dissolved suspended particles) has broad impact in environmental protection and resource utilization via groundwater contamination. The leaching (displacement) of salts and nutrients in soils also has an impact on agricultural production. To predict the contaminant migration in the geological formation more accurately, many analytical solutions for partial differential equations exist but because of the difficulty in obtaining analytical solutions, numerical solutions are more generally used. Numerical solutions are often more difficult to verify, so mathematical model error has to be kept as small as possible. www.intechopen.com

Review of Darcy's law
In 1856, French engineer Henry Darcy was working on a project involving the use of sand to filter the water supply. He performed laboratory experiments to examine the factors that govern the rate of water flowing through the sand. The results of his experiments defined the empirical principle of groundwater flow, in an equation now known as Darcy's law which states that "the saturated flow of water through a column of soil is directly proportional to the head difference and inversely proportional to the length of the column". Darcy's apparatus consisted of a sandfilled column with an inlet and an outlet for water.

Fig. 1. Darcy's law experiment
Two manometers measure the hydraulic head at two points within the column (h1 and h2). The sand is saturated, and a steady flow of water is forced through it at a volumetric rate of Q [L3/T] (Q is sometimes called the volumetric flow rate or the discharge rate). Darcy found that Q was proportional to the head difference dh between the two manometers, inversely proportional to the distance between manometers l, and proportional to the cross sectional area of the sand column (A). This can be mathematically written as: Therefore we can say, Q α A dh/ l Combining these observations and writing an equation in differential form gives Darcy's law for one-dimensional flow: Where Q = volumetric flow rate or the discharge rate(m 3 /s), Cross sectional flow area perpendicular to l(m 2 ), K=hydraulic conductivity (m/s), and d = denotes the change in h over the path l. It can be re-written in differential form as: The minus sign is necessary because head decreases in the direction of flow (i.e., water is always flowing from higher hydraulic head to lower hydraulic head). If there is flow in the positive l direction, Q is positive and dh/dl is negative 2 . Conversely, when flow is in the negative l direction, Q is negative and dh/dl is positive. The constant of proportionality K is the hydraulic conductivity in the l direction, a property of the porous medium and the fluid (water) filling the pores. The common units for hydraulic conductivity are meters/year for regional studies, m/day for local aquifer-scale studies, and cm/sec for laboratory studies. Therefore, in some analysis, we often deviate from the rule of using the SI unit. Another form of the Darcy's law is written for the Darcy flux (or the Darcy Velocity, or, the Specific Discharge) (q) which is the discharge rate per unit cross-sectional area: The Darcy flux q has unit of velocity [L/T] and assumes that flow occurs through the entire cross section of the material without regard to solids and pores. However, Darcy flux is not the actual fluid velocity in the porous media; it is just discharge rate (Q) per unit cross-sectional area. In a Cartesian x, y, z coordinate system, it is commonly expressed as: Same idea applies to estimating conductivities in the y and z directions.

Groundwater velocity / limits of Darcy's law
Application of Darcy's law has both upper and lower limits, it does not hold at very high fluid velocities. Hence, the need to calculate the groundwater velocity because the Specific Discharge (Darcy's flux) of equation 3 which is sometimes called Darcy's velocity is not actually groundwater velocity but discharge rate per unit cross-sectional area (Q/A). Flow of water is limited in pore spaces only, Hence in calculating for the actual seepage velocity of groundwater, the actual cross-sectional area through which the flow is occurring must be accounted for as follows: where V = Groundwater velocity, also known as the seepage velocity, and commonly called average linear velocity (m/s), Q = Volumetric flow rate or the discharge rate(m 3 /s), A= Cross sectional area perpendicular to l (m 2 ),  = porosity.
This is a dimensionless ratio where v = the velocity (m/s), ρ = fluid density (kg/m 3 ), µ = the fluid velocity (kg/m/s), and d = diameter of the pipe (m). Experimental evidence indicates that Darcy's law is valid as long as Re does not exceed a critical value between (1 and 10).This law holds for low velocity and high gradient (Shazrah et al 2008).

Groundwater modelling
Modeling has emerged as a major tool in all branches of science (Igboekwe et al.,2008). Models are conceptual descriptions, tools or devices that represent or describe an approximation of a field situation, real system or natural phenomena. Models are applied to a range of environmental problems mainly for understanding and the interpretation of issues having complex interaction of many variables in the system. They are not exact descriptions of physical systems or processes but are mathematically representing a simplified version of a system. This mathematical calculation is referred to as simulations. Groundwater models are used in calculating rates and direction of groundwater flow in an aquifer. The simulation of groundwater flow needs a proper understanding of the complete hydrogeological characteristics of the site because the applicability, reliability or usefulness of a model depends on how closely the mathematical equations approximate the physical system being modelled. So, the following should be well investigated and understood:  Hydrogeologic framework (surface thickness of aquifers and confining units, boundary conditions that control the rate and direction of groundwater flow).  Hydraulic properties of aquifer( a depiction of the lateral and vertical distribution of hydraulic head i.e initial, steady-state and transient conditions; distribution and www.intechopen.com quantity of groundwater recharge / discharge; and sources or sinks being referred to as stresses and they may be constant or transient). Groundwater modelling entails simulation of aquifer and its response to input/output systems. By mathematically representing a simplified version of a hydrogeological condition, realistic alternative settings can be predicted, tested, and compared. Groundwater models are used to predict/illustrate the groundwater flow and transport processes using mathematical equations based on certain simplifying assumptions. These assumptions normally involve the direction of flow, geometry of the aquifer, the heterogeneity or anisotropy of sediments or bedrock within the aquifer, the contaminant transport mechanisms and chemical reactions. Because of the simplifying assumptions embedded in the mathematical equations and the many uncertainties in the values of data www.intechopen.com required by the model, a model must be viewed as an estimate and not an accurate replication of field settings. Groundwater models, though they represent or approximate a real system are investigation tools useful in many applications. Application of existing groundwater models include water balance (in terms of water quantity), gaining knowledge about the quantitative aspects of the unsaturated zone, simulating of water flow and chemical migration in the saturated zone including rivergroundwater relations, assessing the impact of changes of the groundwater regime on the environment, setting up monitoring networks and groundwater protection zones (Kumar 2002).

Types of models
Models could also be classified based on their typical applications; an illustration is below in table 1.

MODEL APPLICATION Groundwater flow
The problem of water supply is normally described by one equation mainly in terms of hydraulic head. The resultant model providing the solution to this equation is the groundwater flow model.
Water supply, Regional aquifer analysis, Near-well performance, Groundwater/surface water interactions, Dewatering operations.
Solute transport When the problem involves water quality, then an additional equation to the groundwater flow equation is needed to solve for the concentration of the chemical species. The resultant model is the solute transport model. There are other subdivisions of models such as those describing porous media and those of fractured media namely physical scale model, analog model (Hale-Shaw model) etc.

Model concept
Groundwater modelling begins with a conceptual understanding of the physical problem. The next step in modelling is translating the physical system into mathematical terms. In general, the results are the familiar groundwater flow equation and transport equations.

Groundwater flow equation
It is an established hydraulic principle that groundwater moves from areas of higher potential i.e recharge areas (higher elevation or higher pressure/hydraulic head) to areas of lower pressure or elevation. This implies that direction of flow of groundwater ideally follows the topography of the land surface. Cracks, inter-connected pore spaces make a rock material permeable. Some permeable materials may allow fluid to move several metres in a day; while some may move a few centimeters in a century. In the real subsurface, groundwater flows in complex 3D patterns. Darcy's law in three dimensions is analogous to that of one dimension. This is often derived using a representative fixed control volume element (RFCVE) of fixed dimensions usually a cube. of mass) that is 8litres minus 5litres which is 3litres every minute. This implies that at every minute, the water in the RFCVE increases by 3litres in volume. This net flux of mass (change in volume) is known as change in storage (Ss). But on the other hand, the source of the recharge looses to the RFCVE about 8litres of water per minute, thus giving rise to the terms "source" and "sink". Let q be the specific discharge, so that the rate of flow of water through the RFVCE can be expressed in terms of the three components q x, q y, q z. In accordance to the mass balance equation, influx minus outflux is equal to change in storage i.e (q in -q out ) = (Ss). Now considering the flow along x-axis of the RFVCE, the influx through the face, ∆Y∆Z is the same with the specific recharge q x from the x-axis which is q x (x); the outflux is denoted as q This is differentially written as; dqx XYZ dx   = change in flow rate along x axis.
Therefore, change in flow rate along Also, change in flow rate along the y-axis is specified as While that of z-axis is correspondingly given as At this point, the total change in volumetric flow rate is obtained by summing the individual changes in flow rate of the three axes i.e equations (9) The specific change in storage (Ss) is defined as the volume of water released from storage per unit change in head (h) per unit volume of aquifer i.e.

V Ss hXYZ
Then, the rate of change of storage is given by: This could similarly be re-written as is the rate of change in storage, so by bring equation (16)  Equation (16) is the governing equation of a 3D groundwater flow in a confined aquifer because herein, recharge (source) and leakage (sink) is ignored. As earlier on stated, movement of groundwater is based on a hydraulic principle, therefore a source or sink is usually associated within an elemental volume which is expressed as Q, so the general governing equation of a 3D groundwater flow in an unconfined aquifer (which is the basic assumption) in this context is expressed as: where, Kx, Ky, Kz = hydraulic conductivity along the x, y ,z axes which are assumed to be parallel to the major axes of hydraulic conductivity; h = piezometric (hydraulic) head; Q= volumetric flux per unit volume representing source/sink terms; Ss = specific storage coefficient defined as the volume of water released from storage per unit change in head per unit volume of porous material. Equation (19) is the general governing equation of groundwater flow. Groundwater flow in aquifers is often modelled as two-dimensional in the horizontal plane. This is due to the fact that most aquifers have large aspect ratio like a laminated plane sheet of paper, with horizontal dimensions hundreds of times greater than the vertical thickness.
In such a setting, groundwater relatively flows along the horizontal plane, which implies that the z component of the velocity is comparatively small. Therefore, a two-dimensional analysis is carried out in conjuncture with the use of transmissivity, by assuming that h=h(x,y,t) only (h does not vary with z, thus (dh/dz = 0). This simplification of modeling 3D aquifer flow as horizontal two-dimensional flow is called the Dupuit-Forchheimer approximation.
Where, T x and T y is transmissivity in the x and y direction and T = k * b = (m/day)(m) = m 2 /day, while b is the saturated thickness of the aquifer. But a situation, where the aquifer becomes unconfined, b = h, therefore T = k * h , so the equation (20) Where h is the water level function of x and y. Then, the partial differential equation of second order for groundwater flow is given as: where S = storage coefficient; and h = hydraulic head.

Solute transport equation
The transport of solutes in the saturated zone is governed by the advection-dispersion equation which for a porous medium with uniform porosity distribution is formulated as follows: Where, ∂ = delta function meaning change in; x = the dimension (m); t = the time (s) c = concentration of the solute (kg/m 3 ); Rc = reaction rate (concentration of solute) in the source or sink (kg m 3 /s); Dij = dispersion coefficient tensor (m 2 /s); vi = velocity tensor (m/s). An understanding of these equations and their associated boundary and initial conditions is necessary before a modelling problem can be formulated. From equation 23 above, the first term on the right hand side represents advection transport and describes movement of solutes at the average seepage velocity of the flowing groundwater.
The second term represents the change in concentration due to hydrodynamic dispersion. Hydrodynamic dispersion is defined as the sum of mechanical dispersion and molecular diffusion. The third term represents the effects of mixing with a source fluid of different concentration from the groundwater at the point of recharge or injection. Chemical attenuations of inorganic chemicals can occur by sorption/desorption, precipitation/dissolution, oxidation /reduction, etc.
On the other hand, organic chemicals could be absorbed or degraded through microbial processes thus giving rise to a new set of equation such as: Where Y represents all the chemical, geochemical and biological reactions responsible for the transfer of mass between the liquid and solid phase or the conversion of dissolved chemical species from one form to another.
Assuming that the reactions are limited to equilibrium-controlled sorption or exchange and first-order irreversible rate (decay) reactions in an isotropic homogenous porous medium, then equation 24 may be written as: For equilibrium sorption and exchange reaction C t   as well as C is a function of C only.
Hence equilibrium reaction for C and C t   can be substituted into the governing equation to develop a partial differential equation in terms of C alone. The resultant single transport equation is solved for solute concentration, sorbed concentration can be calculated using the equilibrium reaction. The linear sorption reaction considers that the concentration of solute sorbed to the porous medium is directly proportional to the reaction.
where K d = distribution coefficient (M -1 L 3 ) This reaction is assumed to be instantaneous and reversible. The relationship between the adsorbed and dissolved concentrations can be described using three possible isotherms: linear, Langmuir and Freundlich. Thus in terms of the linear isotherm, Isotherms can be incorporated into the transport model using a retardation factor, R [M L -3 ] (Kutílek and Nielsen 1994;Zheng and Bennett 2002): By substituting equation 27 into equation 24, we get Now, substituting equation 28 into equation 29, the resultant equation is as follows:

Boundary conditions
Boundary conditions indicate how an aquifer interacts with the environment outside the model domain. They include things such as heads at surface waters in contact with the aquifer, the location and discharge rate of a pumping well or a leaching irrigation field. For a distinct solution, at least one distinctive boundary condition is specified.
There are four types of boundary conditions which are derived from the most common two: specified head and specified flux conditions. i. Constant Head Boundary: This is a type of specified head boundary condition, in which the head is known and the source of water has a constant water level at the model boundary. This condition is used in modelling an aquifer that is in good interaction with a lake, river or another external aquifer. These are usually where the groundwater is in direct contact with surface water such as a lake or a river and drains interact freely with the aquifer. This condition is also known as the first type boundary condition and it is mathematically referred to as Dirichlet boundary. It is stated mathematically as ii. Constant Flux Boundary: This is a type of specified flux boundary condition also known as the second type of boundary condition, and mathematically known as Neumann's condition or recharge boundaries. Entering or leaving the aquifer is prescribed/constant flux. This boundary condition is used in simulating rainfall or distributed discharge for instance evaporation and also used in specifying known recharge to the aquifer owing to induced recharge or reticulation.
It is stated mathematically as     2 ,N e u m a n n hx hx kx nn is the specified outward normal gradient to the boundary segment 2   .
iii. Mixed Boundary: This type of boundary condition involves some combination of head and flux specification whereby, the rate of flow in and out of the aquifer is a function of the elevation of the stream bed, aquifer head, and leakage between the aquifer and the stream. For instance, leakage through a silty river bed to an underlying aquifer is www.intechopen.com represented by a flux that is proportional to the vertical conductivity of the silt layer and proportional to head difference from the river to the underlying aquifer. It is referred to as the third type of boundary condition and mathematically known as Cauchy Robbins condition. This less common, mixed boundary or induced flux condition is also known as Stream or River Head Dependent Boundary. This boundary condition is used in the modelling of streams in poor connection with the aquifer, upward leakage in artesian aquifers, drains and overlying aquitards. In the figure above, the pumping well is a specified flux condition at the permeable section of the well. The recharge is a specified flux applied at the water table. The low-K bedrock is considered a special specified flux boundary (no-flow boundary). The leaky silt layer below the river is a mixed condition where the flux through the layer is proportional to the difference between the head in the river and the head in the aquifer beneath the silt layer. Where the river is in direct contact with the aquifer, there is a specified head condition.

Mathematical model
Mathematical models are used in simulating the components of the conceptual model and comprise an equation or a set of governing equations representing the processes that occur, for example groundwater flow, solute transport etc. The differential equations are developed from analyzing groundwater flow (and transport) and are known to govern the physics of flow (and transport). The reliability of model predictions depends on how well the model approximates the actual natural situation in the www.intechopen.com field. Certainly, simplifying assumptions are made in order to construct a model, because the field situations are usually too complicated to be simulated exactly. Mathematical models of groundwater flow and solute transport can be solved generally with two broad approaches: Analytical solution of the mathematical equation gives exact solution to the problem, i.e., the unknown variable is solved continuously for every point in space (steady-state flow) and time (transient flow).Analytical models are exact solution to a specified, well simplified groundwater flow or transport equation. Because of the complexity of the 3D groundwater flow and transport equations, the simplicity inherent in analytical model makes it impossible to account for variations in field conditions that occur with time and space. For these problems, (variations in field conditions) such as changes in the rate/direction of groundwater flow, stresses, changes in hydraulic, chemical and complex hydrogeologic boundary conditions, the assumptions to be made to obtain an analytical solution will not be realistic. To solve mathematical models of this nature, we must resort to approximate methods using numerical solution techniques.  Using either FDM or FEM, the partial differential equation is converted to a set of algebraic equations involving unknown values at these particular points. This implies that the partial differential equations representing the balances of considered extensive quantities are replaced by a set of algebraic equations.  Analytical solutions are not available for many problems of practical interest while the numerical methods allow us to solve the governing equation in more than one dimension, for complex boundary conditions and for heterogeneous and anisotropic aquifers, whereas most analytical solutions are restricted to consideration of homogeneous, isotropic aquifers. This implies that the solution is obtained for specified set of numerical values of the various model coefficients.

Finite difference method
The primary vehicle of calculus and other higher mathematics is the function. Its "input value" is its argument usually a point ("P") expressible on a graph. The difference between two points, themselves, is known as their Delta (∆P),as is the difference in their function result, the particular notation being determined by the direction of formation:  Forward difference: ∆F(P) = F (P +∆ P) -F(P);  Central difference:  F(P) = F(P +½∆ P ) -F(P -½∆ P );


Backward difference:  F(P) = F(P) -F(P -∆P). The general preference is the forward orientation, as F(P) is the base, to which differences (i.e "∆P"s) are added to it. Furthermore,  If [∆P] is finite (meaning measurable), then ∆F(P) is known as the finite difference with specific denotations of DP and DF(P);  If [∆P] is infinitesimal (an infinitely small amount -l-usually expressed in standard analysis as a limit;   0 lim  , then ∆F(P) is known as the infinitesimal difference with specific denotations of dP and dF(P) (in calculus graphing,the point is almost exclusively identified as 'x' and F(x) as 'y' The function difference divided by the point difference gives the difference quotient (Newton's quotient): If ∆P is infinitesimal, then the difference quotient is a derivative, otherwise it is a divided difference The approximation of derivatives by finite differences plays a central role in finite difference methods for the numerical solution of differential equations, especially boundary value problems.
An important application of finite differences is in numerical analysis, especially in numerical differential equations, which aim at the numerical solution of ordinary and www.intechopen.com partial differential equations respectively. The idea is to replace the derivatives appearing in the differential equation by finite differences that approximate them. The resulting methods are called finite difference methods. A forward difference is an expression of the form Depending on the application, the spacing h may be variable or constant. A backward difference uses the function values at x and x − h, instead of the values at x + h and x: Finally, the central difference is given by The derivative of a function f at a point x is defined by the limit If h has a fixed (non-zero) value, instead of approaching zero, then the right-hand side is Hence, the forward difference divided by h approximates the derivative when h is small. The error in this approximation can be derived from Taylor's theorem. Assuming that f is continuously differentiable, the error is The same formula holds for the backward difference: However, the central difference yields a more accurate approximation. Its error is proportional to square of the spacing (if f is twice continuously differentiable): The main problem with the central difference method, however, is that oscillating functions can yield zero derivative. If f(nh)=1 for n uneven, and f(nh)=2 for n even, then f'(nh)=0 if it is calculated with the central difference scheme. This is particularly troublesome if the domain of f is discrete.
In an analogous way one can obtain finite difference approximations to higher order derivatives and differential operators. For example, by using the above central difference formula for f'(x + h / 2) and f'(x − h / 2) and applying a central difference formula for the derivative of f' at x, we obtain the central difference approximation of the second derivative of f: More generally, the n th -order forward, backward, and central differences are respectively given by: Note that the central difference will, for odd n, have h multiplied by non-integers. This is often a problem because it amounts to changing the interval of discretization.
Higher-order differences can also be used to construct better approximations. As mentioned above, the first-order difference approximates the first-order derivative up to a term of order h. However, the combination approximates f'(x) up to a term of order h 2 . This can be proven by expanding the above expression in Taylor series, or by using the calculus of finite differences. If necessary, the finite difference can be centered about any point by mixing forward, backward, and central difference.

Basics of finite difference method
The partial differential equation describing the flow and transport processes in groundwater include terms representing derivatives of continuous variables. Finite difference method is based on the approximations of these derivatives (slopes or curves) by discrete linear changes over small discrete intervals of space and time. A situation where the intervals are adequately small, then all of the linear increment will represent a good approximation of the true curvilinear surface. Following an illustration in figure 6a below, Bennett (1976)  These approximations can also be obtained by Taylor series expansion. A certain error is involved in the approximation of the derivatives by finite differences, but this error will generally decrease as ∆x and ∆y are given small values. This error is called 'truncation error' because the replacement of a derivative by a difference quotient is equivalent to a truncated Taylor series. Taylor

Taylor Series and finite difference approximations
The Taylor Series is key to understanding the finite difference method. Considering also the discretization of time, which may be viewed as another dimension and hence represented by another index. Illustrating using a hydrograph in figure 8 below, the head is plotted against time for a transient flow system; n is the index or subscript used to denote the line at which a given head is with respect to time, and it can be approximated as hh tt    . In terms of the heads calculated at specific time increments (or time slope), the slope of the hydrograph at time n can be approximated as follows:   ) at time node t n may be approximated by h t   (Konikow,1996).

The finite-difference grid
The application of a numerical model to the solution of a ground-water problem is a creative process.
There are many different techniques that can be applied to solve the same problem and each modeller has developed preferred ways of approaching a model design and the software package to use in terms of flexibility. However, no software package can be totally flexible. Many software packages are been used by modellers to interactively design a three dimensional finite-difference ground-water flow and contaminant transport models namely GV(Groundwater Vista), MODFLOW, MT3D, and MODPATH. GV model design is generic because it can be used to create data sets for MODFLOW, MT3D, and MODPATH. While each of these specific models has its own data input format, they all have key features in common. The most important features in common are the physical layout of the grid or mesh, the specification of boundary conditions, and the definition of hydraulic properties. A finite-difference model is constructed by dividing the model domain into square or rectangular regions called blocks or cells. Head, drawdown, and concentration are computed at discrete points within the model called nodes. The network of cells and nodes is called the grid or mesh. These terms are used. There are two main types of finitedifference techniques, known as block-centered and mesh-centered. The name of the technique refers to the relationship of the node to the grid lines. Head is computed at the centre of the rectangular cell in the block-centered approach. Conversely, head is computed at the intersection of grid lines (the mesh) in the mesh-centered technique. The figure below illustrates this concept graphically. Fig. 9. The finite-difference grid.
In this figure, note that the dependent variable (head or concentration) is computed at the centre of cells in the block-centered technique but may be offset from the centre in the meshcentered approach.
In each technique, the head and all physical properties are assumed to be constant throughout the cell region surrounding the node. The finite-difference grid is designed by manipulating rows, columns, and layers of cells. A series of cells oriented parallel to the x-direction is called a row. A series of cells along the ydirection is called a column. A horizontal two-dimensional network of cells is called a layer. Cells are designated using the row and column co-ordinates, with the origin in the upper www.intechopen.com left corner of the mesh. That is, the upper left cell is called (row 1, column 1). The upper layer is layer 1and layers increase in number downward. For example, a block-centered finite-difference grid is created in GV by first specifying the number of rows, columns, and layers.
The user also provides the initial row and column widths or spacings. GV then creates a mesh with uniform row and column widths. This is called a regular mesh. While the regular mesh represents the most accurate form of the finite-difference solution (Anderson and Woessner 1992), it is often necessary or desirable to refine the mesh in areas of interest. In this manner, more accuracy is achieved in key areas at the expense of lower accuracy at the edges of the model grid. (Modified after Bear, 1979).
Based on the notation in the figure above, the derivatives at (i, j) for the forward, backward, central and second order of the finite difference formulas of equations (32), (33), (34) and (40) are as follows respectively: at each node i, this is the Implicit method. Implicit method requires much work than Explicit method in solving the set of simultaneous equations, and it is unconditionally stable irrespective of the time step t  .

Crank-Nicholson method
In an Implicit method, the choice of ∆t can be independent of ∆x. Further refinement is possible. Also, the smaller the truncation error, the faster is the convergence of the finite difference equations to the differential equation.
In Implicit and Explicit methods for the time derivative u t   ,the truncation error was of the order of 0(∆t). By replacing u t   with the central difference, the truncation error associated with the time term would be reduced from 0(∆t) to 0{(∆r)²}. This is known as the Crank-Nicholson method.
Modification of equation (65) to handle problems of saturated and unsaturated flow is possible and the extension to more than one space dimension poses no difficulty.

Conclusion
Mathematical models are efficient tools commonly used in studying groundwater systems. Generally, mathematical models are used to simulate (or to predict) the groundwater flow and in some cases the solute and/or heat transport conditions. They may also be used in the evaluation of remediation alternatives. Finite difference methods can be used to solve the solute transport equation, either using backward, forward or central differencing. These methods however can result in artificial oscillation (under or over shooting) or numerical dispersion due to truncation errors of the discretization. Application of backwards finite difference, results in an implicit scheme www.intechopen.com which is always convergent but is computationally expensive, and introduces considerable numerical dispersion (Zheng and Bennett 2002). The use of central differencing (such as Crank Nicholson schemes) in the discretization can cause numerical oscillation in the form of "wiggles" when implicit schemes are used. If an explicit formulation is used instead the solution often is non convergent (Leonard 1979). Numerical oscillation can be minimized by the use of upstream weighting, but this leads to considerable numerical dispersion owing to truncation errors (Zheng and Bennett 2002). Another solution for artificial oscillation is the use of finer grids, with a choice based on the Peclet number: Pe = u *∆x / D where u is the flow velocity [L T -1 ] , Δx is the grid spacing [L] and D is the diffusivity [L 2 T -1 ]. A Pe number < 2 can greatly reduce or eliminate numerical oscillation, but usually the associated computational cost due to excessively fine grids is impractical (Zheng and Bennett 2002). For stability of numerical techniques, the diagonal dominance check is performed, where the coefficient of the left hand side must be equal to the coefficient on the right hand side (Osuwa, J.C 2011).
Because of the associated complexities in numerical modelling, models approximate a real system and therefore are not exact descriptions of physical systems or processes. They're mathematically representing a simplified version of a system, thus expressing a range of possible outcomes reflecting the inherent assumptions. Predictive simulations must be viewed as estimates, dependent upon the quality and uncertainty of the input data. While using models as predictive tools, field monitoring must be incorporated to verify model predictions.
To determine model error, the examination of numerical methods through an analysis of their ability compared with analytical methods is strongly recommended.