King Fahd University of Petroleum & Minerals

Civil Engineering Dept.

Dr. Rashid Allayla

Introductions & Definitions:

Porous Media A medium of interconnected pores that is capable of transmitting fluid. These pores often referred to as “interstices”. Interstices can be fully saturated or partially saturated.

Capillary zone is the zone above the saturated area. It is a zone where pores contain water under negative pressure.

Capillary fringe is the area between the capillary zone and the water table.

Aquifer is any porous formation that can transmit water at economic (usable) quantity. It comes from the Latin words aqua (means water) and ferro (to bear).

Aquiclude is a saturated formation that is relatively impermeable and not capable of transmitting water under ordinary hydraulic conditions. Clay belongs to this type.

Aquifuge is a formation that neither contains water nor transmits it. Solid granite belongs to this type.

Aquitard is a saturated formation with insufficient permeability to allow completion of production wells. If Aquitard is large enough, it can be an important groundwater storage zone.

Sandstone is sedimentary rock formations that are made of segmented sand and gravel. Non-compacted sandstone has porosity in the range of 30 to 50%, which is highly permeable.

Carbonate rocks are formations that contain gypsum and rock salts. The deeper the formation the more crystallized the minerals and the less permeable. The porosity range of carbonate rock is 20 to 30% but fractured rock could be much higher.

Classification of Aquifers:

A confined aquifer also known as artesian aquifer is a formation containing transmittable quantities water that is under pressure greater than atmospheric pressure. Confined aquifer is bounded from above and from below by impermeable formations. The pressure at this type of aquifers is greater than atmospheric. Unconfined aquifer, also known as water-table aquifer is an aquifer with free water surface that is open to atmosphere. The upper surface of the zone of saturation is the water table. It is the upper portion of the saturated zone and is open to atmosphere. Unconfined aquifers that contain Aquitard and and/or Aquiclude may contain additional water tables.

A piezometric surface is the contour of water elevation of wells tapping confined aquifers. The contours provide indications of the direction of groundwater flow in the aquifers. If the piezometric surface falls below the top formation of confined aquifers (as shown in the left figure below), the aquifer at that point is a water-table aquifer and the surface becomes water-table surface. It is easy to confuse water-table surface with piezometric surface when both confined and unconfined aquifers exist on to of each other. But , in generals the two surfaces do not coincide. A piezometer is a device, which indicates the water pressure head at a point in the confined aquifer. It consists of casing that is open in both sides and fits tightly against the geological formation making up the aquifer. The height to which water rises in the piezometer is the water pressure head.

Averaging Process in Porous Media:

If for a given point P in a porous medium, the volume that contains fluid is V_v and the total volume of the porous medium is V, the porosity is the ratio between the volume of voids and the total volume of the porous medium V. For a large value of V, porosity φ is constant unless the porous media is inhomogeneous. Reducing V beyond a certain value (say U₀), φ starts to fluctuate. This happens as V approaches the dimensions of a single pore. Reducing V even further, φ will become either one or zero depending on whether the point is inside a void space or a solid matrix. The point (P) at which the ratio φ become meaningful is known as the representative elementary volume (REV). At this point,

∂φ(P) / ∂V₀ = 0 for 0 < φ > 1

REV

Porosity: It the ratio of the interconnected voids within the soil to the total volume of the soil. If V_V designates the volume of “interconnected” voids, V_S is the bulk volume and V_T is the total volume of the soil media, the porosity θ is defined as:

θ = V_V / V_T = (V_T-V_S) / V_T = 1 – V_S/V_T

The above definition is also referred to as the effective porosity.

EXAMPLE:

A Soil sample was taken from a cavity in the field. The sample weighs 200 grams when it was at field moisture. The depth of the cavity is 5 cm and volume is 130 cm³. After drying the sample, it weighed 180 grams. The entire weight of the sample under water is 111 grams. Find: 1) Dry bulk density, 2) Wet bulk density, 3) Volumetric water content, 4) Total porosity and 5) Degree of Saturation.

1) Dry bulk density: M_d / V = 180/130 = 138 gm / Cm³

2) Wet bulk density: (M_d + M_w) / V = [180-(200-180)] / 130 = 1.54 gm / Cm³

3) Volumetric water content = V_w / V = (200-180) / 130

Volume of solid = Volume of displaced liquid

Actual weight – weight of water = buoyancy = 180 – 111 = 69 gm, therefore V_s= 69 cm³

So the particle density = M_s / V 180 / 69 = 2.61

4) Total porosity = V_p / V = (130-69) / 130 = 0.47

5) Degree of saturation = θ / θ_s = (V_w / V_b) / (V_p / V_b ) = 0.15 / 0.47 = 0.32

Measurement of Porosity

1. Direct Method: This method is accomplished by crushing the porous media sample to determine the volume of the solid. This method is crude since it removes non-interconnected pores, which, in groundwater analysis, are considered part of the solid volume because they do not contribute to the flow.

2. Mercury Injection Method: This is accomplished by placing a porous medium sample in mercury inside vacuum chamber. Mercury is forced inside the pores by increasing the pressure of air in the chamber. The method is not suitable for porous medium with low permeability, as it requires very high pressure to force Mercury into the small pores.

Porous Medium

Typical Values of Porosity:

Material	Porosity
Clay	0.45-0.55
Silt	0.35-0.50
Sand	0.25-0.40
Gravel	0.25-0.40
Limestone	0.01-0.20
Shale	0.0-0.10
Fractured Rock	0.0-0.10
Dense, Solid Rock	<0.01

Distribution of Pressure:

The force that balances the gravity force and holds water in equilibrium is,

∂P / ∂ z = - γ

Where P is the pressure and γ is the unit weight of water. Integrating,

P = - γz + C

Selecting water table as the datum, pressure above water table is negative and pressure below is positive. One can ask how does water stay in voids and not drain. The answer is because of Surface Tension.

Capillary Action:

Capillary action is the result of surface tension acting between water particles and solid particles. It occurs when adhesion to the walls of a small tube (or pores in porous media) is stronger than cohesive forces between water molecules. Adhesion of water to the walls of a vessel will cause an upward force on the liquid at the edges and, as a result, an upward movement of water takes place. The height of the water column produced depends on the magnitude of surface tension σ between fluid and solid phases. When the system reaches equilibrium, the force of the interfacial tension is equal and opposite in direction to the weight of water.

Pressure = - γ h

The downward force = (-γh) (πr²)

The Upward force = (Surface tension) (circumference of the tube)

F_Upward = F_Downward

(-γh) (πr²) = σ Cos α (2πr)

Then, for small α:

h = (2σ) / γ r

Introducing the term capillary pressure (P_c) which is defined as the negative pressure that water particles are under in the partially saturated zone located above water table. The capillary pressure head, P_c/γg is equivalent to the distance h above the water table. The above equation indicates that the curvature of the interface changes with P_c (or h). Since P_c increases with elevation, the radius of curvature of the interface decreases with elevation. So, as P_c is increased, the less saturated soil is.

Fundamentals of Flow Equations:

Reynolds’s Transport Theorem:

Flow rate in – flow rate out =  V₂A₂ – V₁A₁ = V₂ . A_{2 –} V₁ . A₁ = Σ V . A

The mass rate of flow out of control volume is,

m = Σ V. A

The rate of flow of property B (property B could be mass, energy or momentum per unit time) out of the control volume is the product of mass rate m with intensive property β defined as B per unit mass)

B = Σ β ρV . A

Where, ρVA is the mass rate. To verify the above

Mass = (mass/mass). (mass/L³⁾. (L/time). (L² )

If the velocity varies, it is necessary to integrate the velocity over the area A.

B = β ρ ò V. dA      

Reynolds Transport Theorem States that,

[TIME RATE OF CHANGE OF N (SAY MASS) OF A SYSTEM] = [TIME RATE OF CHANGE OF THE CONTENT OF CONTROL VOLUME] + [NET RATE OF N THROUGH THE CONTROL SURFACE]

Mathematically, it states that,

D / Dt ò_S ρdV = ∂/∂t ʃ_CV (ρ)dV + ʃ_{CS (}ρV.n) dA

D/Dt is the material derivative. Also If V. n > 0, the property is carried into the control volume. If V. n = 0, no transport of N is carried out of CV.

For steady Conditions, the equation become,

D/Dt ʃ_CV ρdV = ∂/∂t ʃ_CV ρdV + ʃ_CS ρV.n dA = 0

= 0 = 0

which states that the mass rate of flow out of a region in space minus the mass rate of flow out of the region is equal to zero if the time rate of change in mass is unchanged. The equation becomes, (ρ VA) _out – (ρ V A) _in = 0

The mass continuity equation for incompressible fluid is,

∂q_x / ∂x + ∂q_y / ∂ y + ∂q_z / ∂z = 0

EXAMPLE:

A rectangular cross sectional bathtub is being filled with water from a faucet. The steady flow rate is .01 cubic meter per minute. Estimate the time rate of change of depth of water ∂h / ∂t in cm/min at any given instant

Faucet

h

The continuity, ∂/∂t ∫_CV ρdV + ∫_CS ρV.n dA

Here: V.n dA is equal to – ρQ

The first term: ∫_CV ρdV where ∫dV is = (h) (2m)(5 m)=10 h then(10) (ρ∂h / ∂t) – ρ Q = 0 Then,

∂h / ∂t = Q/10 = 0.01 cu meter per min/10 cu meter = 0.01 cm/min

For a compressible medium (such as baloons), the continuity equation is:

(∂q_x / ∂x + ∂q_y / ∂ y + ∂q_z / ∂z) (∆x ∆y ∆z) = Change in storage

In Groundwater, the change of storage is represented by specific storage S_s which is defined as the volume of water released from storage per unit change in head per unit volume of aquifer:

Ss = - ∆V /(∆x ∆y ∆z)

The convention here is that ∆V is positive when ∆h is negative and vice versa. The rate of change in storage is ∆V/∆t, then:

∆V/∆t = - Ss (∆h / ∆t) (∆x ∆y ∆z)

Combining with the above and dividing through by (∆x ∆y ∆z), the ground water balance equation become:

(∂q_x / ∂x + ∂q_y / ∂ y + ∂q_z / ∂z) = - S_s (∂h / ∂t) S_s measured in 1/L

The above equation has theoretical meaning but very little practical application because q cannot be measured accurately and it must be related to a specific discharge (measurable term) through the use of Darcy’s law.

Unit Cross Sectional

Area

S_s

Elements of Water Retention Curve:

Capillary Fringe: It is a distance immediately above water table that is saturated but under negative pressure. The absolute value of this pressure is not large enough to displace the water.

Specific Yield (S_Y): It is the fraction of the saturated bulk volume of water, which will drain by gravity when the water table drops. Specific yield is less than porosity because some water adsorbed too strongly to the soil particles to drain.

Specific Retention (S_r) of a rock or soil is the ratio of the volume of water a porous medium can retain against gravity to the total volume of the porous medium.

S_Y = φ - S_r

Note that measuring specific yield in the field; the value is always less that the value measured from the formula above because, in the field, there is no instantaneous release of water.

Apparent Specific Yield S_ya:

This parameter is the ratio of the volume of water drained from the saturated aquifer to the resulting change in the volume of the aquifer. Because drainage of saturated pores is not instantaneous, the measured apparent specific yield does not reach its upper value of S_y. When the time lag of complete drainage is taken into account, it takes some time for drainage to be complete and the specific yield become time dependent and the value taken at any time is the apparent specific yield. The upper limit of the apparent specific yield can be reached only when the aquifer is homogeneous, drainage is complete and the fall of water table is much greater than the height of capillary fringe.

Position 1

Position 2

Volumetric Water Content (θ) is the fraction of porous media that contains water.

Water-Retention Curve is capillary pressure-desaturation curve, which display the ability of porous medium to retain water above water table when water is under negative pressure.

Figures below show typical water retention curves. At one region in the curve (immediately above the water table), the negative pressure is not sufficiently large to desaturate the porous medium (area known as capillary fringe region). Significant desaturation occurs only when the difference of pressure and air is high enough to displace water (that is, draining it to the water table). The second figure indicates that, the finer the soil is, the larger the fringe is. Note that there will be a point further along the curve where drainage cannot take place with further increase of capillary pressure head (area is known as water retention region). In agricultural engineering, this region is called the field capacity of soil. As the second figure below shows, fine textured materials contain significantly more pores and tend to retain more water.

S_r

θ

Medium Compressibility:

There are two types of compressibility within an aquifer system:

Internal Stress (hydrostatic pressure) P: Is the pressure of water below saturated zone.

External Stress σ: Is the stress due to overburden pressure. σ = σ^’ + dP where σ^’ is the effective stress or the stress between solid and fluid. A negative stress component indicates tension and positive component indicates compression. In general, stress in aquifer is a symmetrical tenser,

σ_ij = σ^’_ij + P δ_ij where the dirac function indicates that the pressure is an isotropic quantity. Any incremental change in σ produces an immediate change in P because σ_’ responds slowly to changes in σ. Then

dσ = d σ_’ + dP

In the above equation, withdrawing water from a confined aquifers, would lower pore-water pressure. Since the total overburden pressure remains unaffected, any change in pore water pressure is accompanied by a corresponding increase in effective stress. This results in compaction of the porous skeleton and a reduction in the pore volume of the medium. Lateral deformation of the aquifer is generally negligible with respect to vertical deformation and, therefore, will be ignored in future considerations.

P

_σ’

Defining the coefficient of compressibility of water as,

β = - 1/V_w ∂V_w / ∂P = 1/ρ ∂ρ/∂P

The minus sign indicates that as the pressure increases, the volume decreases. Since β is independent of pressure, the volume of water become,

V_w = V_w0 exp [-β(P – P₀)]

Where subscript 0 indicate values at reference pressure P₀.

Similarly, the coefficient of compressibility, α, which is a measure of the elastic property of the soil matrix, is defined as,

α = -1/V_b ∂V_b / ∂/σ^’ where V_b is the volume of the soil matrix.

Considering only vertical deformation, the total volume of porous soil matrix V_b is V_s/(1-φ) and since ∂V_s/∂σ^’ = 0, α becomes,

α = 1/(1-φ) ∂φ/∂P which relates α to the change of porosity φ.

Specific Storage, S_s and Storage coefficient, S:

It is important now to find a way to calculate the amount of water released from the aquifer when the pressure head is changed. When we reduce water pressure by withdrawing water from a confined aquifer, water is released by to two mechanisms. The first release is due to the compressibility of water, represented by β and the second release is due to the compressibility soil matrix, represented by α. Writing the expression for total mass in a saturated element of aquifer as M = ρ φ V_b and assuming aquifer deforms only in Z direction, also, since V_b is ∆X∆Y∆Z, the change of mass is,

dM = [ρd(φ∆Z) + φ∆Zdρ] ∆X∆Y

The first quantity is the contribution per unit area due to change in pore volume at constant ρ. The second is due to change of water density at constant pore volume.

Confined aquifer

From the previous considerations, we can now define the specific storage S_s of a confined aquifer as the amount of water released from the aquifer per unit Volume per unit decline of piezometric head or,

S_s = ∆V_w /V ∆h

Where h is the decline of head. In terms of β and α, the Specific storage become,

S_s = ρgφ (α + β)

Another parameter used for confined aquifer analysis is the storage coefficient S. It is defined as the product of S_s with the aquifer thickness b. This coefficient is dimensionless and defined as the volume of water released from a column of unit area and height b per unit decline of pressure head.

The storage coefficient for unconfined (phreatic) aquifers is defined as the volume of water released per unit area of aquifer per unit decline of water table. Unlike the release mechanism, defining confined aquifers (which are due to compressibility of the medium and water), water drained from the volume of pore space between initial and final positions of water table is due to gravity. The values of α and β play a very little role in drainage for this type of aquifers. The storage coefficient of phreatic aquifer is equivalent to specific yield discussed earlier. Note that because there is residual fluid retained in the pore space (named earlier as specific retention) specific yield cannot be confused with porosity, as the drainage is never a complete one. Therefore,

S_y = φ - S_r

For this reason, S_y is referred to as the effective porosity or drainable porosity

Typical values of S:

Confined aquifers: S ranges from10^-4 to 10^-6 of which 40% due to expansion of water and 60% due to compression of soil

Unconfined aquifers: S_y may be 20% to 30% (From Bear)

What affects specific yield?

Specific yield is a measure of how much water can drain away from the rock under gravity versus how much water the rock actually holds. Since surface retention is proportional to the water-holding capacity of soil particles, grain size plays an important role in determining the value of specific yield. The smaller the particle, the larger of surface area, the less the specific yield. This is the reason why pumping from sandy aquifer yields more water than clayey aquifer.

EXAMPLE: (From McWhorter 1984)

The average volume of a confined aquifer per km² is 3 x 10⁷. The storage coefficient of the aquifer at a location where the thickness b = 50 m is 0.0034. Estimate the volume of water recovered per km² by reducing the pressure head by 25 meters.

From the definition of specific storage as the volume of water released from storage per unit volume of aquifer per unit decline of pressure head. The specific storage is,

S_s = S / b = 0.0034 / 50 = 6.8 x 10^-5 per m

The volume of water recovered per km² = (S_s) (Volume)(head drop) = (6.8 x 10^{-5) (}3 x 10⁷)(25)

= 5 x 10⁴ m³

Groundwater Motion:

Darcy’s Empirical Equation:

The French engineer Henry Darcy introduced Darcy’s Law in 1856 when he was investigating flow under foundations in the French city of Dijon.

The law is a generalized relationship between the flow of fluid and the change in head in porous media under saturated conditions. Darcy concluded that Q (measured in volume of water per unit time) is directly proportional to the cross-sectional area of the porous medium transmitting water, the difference in head between the entrance and the exit and inversely proportional to the length separating the points. This law is valid for any Newtonian fluid and, with adjustment, could be used to accommodate unsaturated flow conditions as well.

Porous Ceramic

For a one dimensional flow perpendicular to the cross sectional area, Darcy found that,

Q = -A K ∆h / L

Where the minus sign because ∆h = h₂ – h₁ is negative and Q must be positive. Q is the volumetric flow rate (m³/s or ft³/s), A is the flow area perpendicular to the flow direction (m² or ft²), K is the hydraulic conductivity of the porous medium (m/s or ft/s), L is the flow path length (m or ft), h is the hydraulic head (m or ft), and ∆h is the change in h over the path L. The hydraulic head at a specific point, h is the sum of the pressure head and the elevation, or

h = (p/γ + z)

Where, p is the water pressure (N/m², lb/ft²), γ is the specific weight of water (lb/ft³), g is the acceleration of gravity (m/s² or ft/s²), and z is the elevation (m or ft).Note that the hydraulic head is the height that water would raise in a piezometer. h is simply the difference in height of water in a piezometer placed at the inlet and the outlet.

Substituting the last formula into the first yield,

Q = AK ∆ [(P/γ) + z] / L

This is the Darcy’s formula for one-dimensional flow. Hydraulic conductivity K

Darcy’s Formula in Three-dimensional Form:

q = -K Del h where Del h is ∂h/∂x i + ∂h/∂y j + ∂h/∂z k, q is the vector velocity discharge equal to Q/A. Or,

Q = ʃ_A n.q dA

Introducing the term intrinsic permeability, which relates to hydraulic conductivity as?

K = k γ / μ

From the above definition, k has the dimension L². Intrinsic permeability pertains to the relative ease with which a porous medium can transmit a liquid under a hydraulic or gradient. It is a property of the porous medium and is independent of the nature of the liquid or the potential field. Note that the hydraulic conductivity K is dependent on both soil medium (represented by k) and fluid properties (represented by μ & γ). Various formulas relate k to properties of porous medium. One of them is,

k = C d² k in cm² and d in cm

C is a constant that varies from 45 for clay sand to 140 for pure sand.

Units:

The standard unit used by hydrologists for hydraulic conductivity, K is meter per day. In laboratory, the standard unit is in gallons per day through area of a porous medium measured in ft². In the field, hydraulic conductivity is measured as the discharge of water through a cross-area of an aquifer one foot thick and one mile wide under a hydraulic gradient of 1 ft/mile (Bear 1979).

1 US gal per day ft² = 4.72 x 10^-5 cm / sec = 4.08 x 10^-2 m/d

The standard unit of permeability k is cm². Other units used is darcy which is defined as,

1 darcy = 1 cm³/sec/cm² x 1 centipoise per one atmosphere per cm

EXAMPLE:

The hydraulic conductivity for a given porous medium is 4.8 x 10^-4 cm/s for water with a density of 1 gm/cm³ and viscosity of 1 centi-poise. Calculate the intrinsic permeability and the hydraulic conductivity for oil with ρ = 0.73 gm / cm3 and μ = 1.8 centi-poise.

K = μ / ρg K = [(0.01 dynes-s/cm²) (4.8 x 10^-4 cm/s) / 980 dynes / cm³ = 4.9 x 10^-9 cm²

K_oil = (4.9 x 10^-9 cm²⁾ (0.73 gm / cm3) (980) / 0.018) = 1.95 x 10^-4 cm / s

Laboratory Measurement of Hydraulic Conductivity:

Measurement of hydraulic conductivity through field aquifer tests provides more reliable results than in laboratory tests because, field tests involve large magnitudes of porous medium compared to small laboratory samples. Field tests will be discussed when aquifer testing method is covered. In laboratories, hydraulic conductivity is measured using Permeameters. Two commonly used permeameters are shown below; one is constant head where the difference between water levels reflects the head loss between inlet and outlet of the soil sample. In the falling-head permeameters, the rate of discharge through the sample decreases with time as the driving head decreases. The equations used for the two measurements are,

K = QL / ∆h A for the constant head permeameters, and,

K = aL / At ln ∆h₀ / ∆h (t) for the falling head permeameters.

Darcy’s law shown above is limited to one –dimensional flow of homogeneous incompressible flow. Introducing the concept of specific discharge q, Darcy’s law becomes,

q = K Δh / L

In the above, q is also referred to as the Darcy flux. It is fictitious form of velocity because the equation assumes that the discharge occurs throughout the cross sectional area of soil in spite of the fact that solid particles constitute a major portion of the cross sectional area. The portion of the area available to flow is equal to φA, the average “real” velocity is, therefore is,

v = q / φ

Specific discharge

Intrinsic velocity

The energy loss Δh is due to the friction between the moving water and the walls of the solid. The hydraulic gradient is simply the difference between the head at inlet and head at the outlet divided by the distance between the inlet and the outlet. Note that the flow prescribed by Darcy’s law,

1) Takes place from higher head to lower head and not necessarily from higher pressure to lower one.

2) When the medium is non-homogeneous, Darcy’s law is still valid since K(x,y,z) lies outside the gradient of head (i.e. q = K (x,y,z) Ñ h).

3) When the medium is anisotropic, K can no longer be considered scalar but tensor causing the vector q and Ñh to be non-collinear (i.e. Not in the same direction)

4) Darcy’s law specifies linear relationship between velocity and hydraulic head. This relationship is valid only for small Reynolds number. At high Reynolds number, viscous forces do not govern the flow and the hydraulic gradient will have higher order terms.

5) Darcy’s law is based on the concept of non-slip conditions at the solid boundary. When fluid is under low pressure or when gas is used as fluid through porous medium, the fluid molecules are not always in contact with the solid particles causing a finite velocity at the solid boundary and Darcy’s equation could not be applied.

q

Non-Homogeneity (Heterogeneity) and Anisotropy:

A medium is non-homogeneous when elements such as hydraulic conductivity vary with space. Rarely is an aquifer actually homogeneous, and due to the difficulties in solving non-homogeneous aquifer system, flow nets are often employed to transform such system to before employing such solutions. A medium is an anisotropic if the hydraulic conductivity is directional. In non-homogeneous aquifers, Darcy’s equation is still valid because the value of K(x,y,z) is still a scalar and lies outside the head gradient. Non-homogeneity in aquifers is often the result of stratification of the aquifer. The individual stratum may be homogeneous but the values of K may vary between one layer to the other making the whole system non-homogeneous

Equivalent Hydraulic Conductivity in Stratified Aquifers:

Consider the stratified aquifer shown. The average value of K can be calculated by treating the stratified aquifer as a single homogeneous medium. Since the head drop across each layer is the same, the discharge rate through the entire stratified confined aquifer is,

Q = B Ќ (h₂ –h₁) / L

Where B = b₁ + b₂ + …. + b_n n is the number of layers, b is the thickness of individual layer and Ќ is the equivalent hydraulic conductivity of the whole system. The equivalent hydraulic conductivity become,

Ќ = Σ b_i K_i / Σ b_i i = 1 to n

In a similar way, if we consider flow through vertical strata in series, the discharge is the same across each strata and the total difference of head is the sum of ∆h between each strata. Since,

∆h_L= ∆h₁ +∆h₂ + …… + ∆h_n

Then,

QL / B Ќ = QL₁ /b K₁ + ….. QL_n /b K_n

And,

Ќ = L / [ L₁ / K₁ + …. + L_n / K_n]

Where,

L = L₁ + L₂ + …. + L_n

b_n

L_1, L₂ … L_n

Anisotropic Aquifers:

When an aquifer is stratified (or bedded) microscopically, and the amount of stratification is smaller than the represented volume element of soil, stratification does not cause the aquifer to be non-homogeneous but it becomes anisotropic. The hydraulic conductivity in the direction parallel to the bedding is generally greater than in the vertical direction (K_x > K_y). Under these conditions, the aquifer’s hydraulic conductivity is said to have directional property. The primary cause of anisotropy is the orientation of clay minerals in sedimentary rocks and can also be the results of fracture of material composing the aquifer. The hydraulic conductivity under these conditions is not scalar quantity and is not collinear with the velocity vector. In general, the hydraulic conductivity is considered as second rank tensor, and if coordinate system is set up in such a way that the coordinate system coincide with the principal directions of anisotropy. Fortunately, in the field conditions, the direction of bedded formations coincides with the principal directions of the hydraulic conductivity.

EXAMPLE: (From McWhorter 77)

One-dimensional flow between two parallel channels occurs in a homogeneous aquifer of thickness 4 m. The difference in head is 1.3 m; the discharge is 1.82 x 10^-5 m³ / S per meter of length of channel. The channels are 10 m apart. A layer of sediment is ultimately deposited at the inflow face with thickness of 4 cm and hydraulic conductivity of 1.4 x 10^-5 cm/S. Calculate the discharge rate per Km of the channels following the deposition of the sediment.

Q = Q/b = 1.82 x 10^-5 m³ / S / 4 m = 4.55 x 10^-6 m/S

The hydraulic conductivity of the aquifer is

K = qL/∆h = (4.55 x 10^-6 m/S) (10 m) / 1.3 m = 3.5 x 10^-3 cm / S

Ќ = (10 + 0.04) / [(0.04 / 1.4 x 10^-5) + (10 / 3.5 x 10^-3 )] = 1.76 x 10^-3 cm/S

The discharge rate per Km between the channels after deposition of the sediment is,

Q = Ќ b (h₂ – h₁) / L = (1.76 x 10^-5) (4m) (1.3m) / 10.04 m x 1000 m/km = 9.1 x 10^-3 m³/S-Km

EXAMPLE:

A permeameters consisting of two layers. The top layer is 10cm in thickness and the bottom layer is 5 cm. The hydraulic conductivity of the top layer is 0.0515 cm/s and the bottom is 0.0022 cm/s. the top and the bottom is subjected to a differential head of 25 cm. Find the intrinsic velocity across the soil. If the top soil has φ = .2

Q = K Δh /ΔL A

Δh = 25cm

Q₁ = Q₂

(0.0515) (Δh₁ / 10) A₁ = (0.0022) (Δh₂ / 5) A₂

But A₁ = A₂ Then Δh₂ = 11.705 Δh₁

Since Δh₂ + Δh₁ = 25

Then Δh₁ = 25 / 12.705 = 1.97 cm

10 cm 5 cm

q = (0.0515) (1.97/10) = 0.0101 cm / s

v = q / φ = 0.0101 / 0.2 = 0.202 cm / s

EXAMPLE:

A tube has a cross sectional area a = 4 cm² contains soil of K = 1000 cm/day and length L – 40 cm. The upstream piezometer has cross section of 8 cm² and piezometric head of 20 cm. The downstream piezometer has cross sectional area of 16 cm² and piezometric head of 30 cm. At time zero, the flow across the tube is allowed to take place. Find the time it takes for upstream and downstream piezometric heads to be at the same level.

L = 40 cm

a = 4 cm²

30 cm

From continuity: A₁ y = A₂ (10 – y)

8 y + 16 y = 160 Then y = 160 / 24 = 6.66 cm

From Darcy’s equation Q = K a (Δh / L)

But Q = amount of volume of water enters and exits the soil in tube per time

(A₁) (y) is the volume

Since Q= Ka Δh / L then

Q = [(8 cm²) (6.66 cm) / t] = (1000 cm/s)(4 cm²) (30 – 20) / (40 cm) therefore t = 4603.4 sec =

Differential Equation of Groundwater Flow:

The equation describing groundwater flow is derived by combining the flux equation described by Darcy’s equation with the equation of mass balance. Consider an elemental volume ∆V (∆x∆y∆z). The net inward flux of water through the elemental volume ∆V must equal to the rate of accumulation within the elemental volume.

Elemental volume

Outflow rate – inflow rate in x-direction = ∂/∂x (ρQ)∆x

The volume discharge Q is the product of Darcy velocity with the cross sectional area normal to flow Q_x = q_x ∆y∆z, Q_x = q_y ∆x∆z, Q_z = q_z ∆y∆x

So the time rate of change of mass,

∂M/∂t = - [∂(ρq_x)/∂x + ∂(ρq_y)/∂y + ∂(ρq_z)/∂z] ∆x∆y∆z

Since ∂(ρq) /∂x is ρ∂q/∂x + q∂ρ/∂x .It was shown before that dρ = -ρ β dP substituting for dρ,

∂M/∂t = -ρ [∂(q_x/∂x)+(q_y/∂y)+∂(q_z/∂z)] ∆x∆y∆z – [q_xρβ ∂P/∂x-q_xρβ ∂P/∂-q_xρβ ∂P/∂x] ∆x∆y∆z

The last three terms are very small since β is very small, the mass balance equation become,

∂M/∂t = - ρ [∂(q_x/∂x)+(q_y/∂y)+∂(q_z/∂z)] ∆x∆y∆z

Introducing the Darcy’s equation in a non homogeneous anisotropic aquifer in which the principal directions coincide with the directions x,y and z, the equation become,

(∂M/∂t) / ρ∆x∆y∆z = ∂/∂x (K_x∂h/∂x) + ∂/∂y(K_y∂h/∂y) + ∂/∂z(K_z∂h/∂z)

From previous analysis, it can be shown that,

∂M/∂t / ρ∆x∆y∆z = S_s ∂h/∂t then

∂/∂x(K_x∂h/∂x) + ∂/∂y(K_y∂h/∂y) + ∂/∂z(K_z∂h/∂z) = S_s ∂h /∂t

The above equation is linear partial differential equation describing the head distribution with respect to time and space of non-homogeneous, anisotropic confined aquifers. When the confined aquifer is homogeneous but anisotropic, the equation become,

S_s ∂h /∂t = K_x ∂²h / ∂x² + K_y ∂²h / ∂y² + K_z ∂²h / ∂z²

And if, K_x = K_y = K_z = K, The equation become,

∂²h / ∂x² + ∂²h / ∂y² + ∂²h / ∂z² = S_s / K ∂h/∂t

The above is the equation describing flow in homogeneous, isotropic confined aquifer. For aquifer of constant thickness, the z term drops yielding,

∂²h / ∂x² + ∂²h / ∂y² = S/bK ∂h/∂t

The product bK is known as the transmissibility (or transmissivity) T of confined aquifer. The term has the dimension of L²/t. Introducing drawdown s at a point in an aquifer where,

S = h₀ – h h₀ is a reference value of piezometric head. The equation becomes,

∂²s / ∂x² + ∂²s / ∂y² = (S/T) ∂s/∂t

If the replenishment rate of the confined aquifer is equal to the outflow rate, the change of storage is zero and the equation takes the form of Laplace equation, which describes the flow in confined aquifer under steady state conditions. The Laplace equation in three dimensions is,

∂²s / ∂x² + ∂²s / ∂y² + ∂²s / ∂z² = 0 or Ѳ s = 0

There are numerous solutions of this equation in math and groundwater literature. Some of the method will be discussed later.

Groundwater Equation for Unconfined Aquifers:

Recall that the volume of water derived from storage of confined aquifers is due to two main processes: expansion of water (as defined by β) and compaction of the soil matrix (as defined by α). On the other hand, water derived from unconfined aquifers is due to three processes: gravity induced drainage of pores, expansion of water and compaction of the soil matrix. The contributions from the latter two processes are minor compared to the first. Thus, we will assume that the gravity drainage will be the only process of water release from unconfined aquifer. The amount of water released is measured by calculating the volume of cone of depression surrounding the pumping well multiplied by the apparent specific yield of the aquifer.

Dupuit Assumption:

The equation describing flow of water in unconfined aquifer is described by the same equation in the bottom of page 22. For unconfined aquifers, S_s is very small and the equation become,

∂/∂x(K_x∂h/∂x) + ∂/∂y(K_y∂h/∂y) + ∂/∂z(K_z∂h/∂z) = 0

Note that the right hand side of the equation is set to zero because S_s in unconfined aquifer is negligible not because the flow is steady. The solution of the above equation yields the location of the water table h(x,y,z,t) which is function of position and time. However, the location of the water table is required (a priori) before one can proceed to solve the equation. The difficulties associated with solution of this equation lead us to seek some approximations. This is known as the Dupuit-Forchheimer approximation.

h = f (x , z)

h = f (x)

The assumption states that, if the slope of the phreatic surface is small and the depth of water is large, the x-component of the specific discharge at the water table is not significantly different from that at the bottom of the aquifer. this implies that the equipotential lines are vertical and streamlines are horizontal which also implies that the flow is horizontal. Since,

q_s = –K dh/ds = K sin θ

Dupuit assumption is equivalent to replacing sin θ (dh/ds) by tan θ = dh/dx. The discharge under this condition become,

Q = -K (x) h dh/dx = K/2L (h₀² – h_L²) measured per unit width of aquifer.

Note that since flow near seepage face (at the outlet) is vertical, Dupuit assumption does not apply at this location.

The value h in the equation represents both the thickness of the flow as well as the piezometric head at the water table. Dupuit assumption of horizontal flow permits the use of mass balance through control volume that extends from the water table to the aquifer bottom. Since water compressibility is not important in unconfined aquifers, we will employ volume balance instead of mass balance to derive the equation describing flow in unconfined aquifers.

Groundwater Equation for Unconfined Aquifer Based on Dupuit Assumption (Boussinesq Equation):

Consider an elemental volume in an unconfined aquifer extending from the water table down to the impermeable boundary as shown. The material balance states that,

∂Q_x/∂x ∆x +∂Q_y/∂y ∆y = net volume of outflow rate (=∂V_w/∂t)

Substituting Darcy’s formula for Q, and if the thickness of the aquifer is unity,

- ∂/∂x (Kh∆y ∂h/∂x)∆x - ∂/∂y (Kh∆x ∂h/∂y)∆y = net volume of outflow rate

Where h∆x is the area perpendicular to flow direction in the x axis and h∆y for the y direction as shown

∆x

x

∆y

Q_x

Q_y

WT at t = 0

WT at t = t

Net outflow rate/∆x∆y = -∂/∂x (K h ∂h/∂x) - ∂/∂y (K h ∂h/∂y)

From the definition of apparent specific yield “The volume of water drained per horizontal area of aquifer per unit decline of head”, then,

∂V_w / ∂t = S_ya ∂h/∂t ∆x∆y

Combining the two equations yield,

∂/∂x (K h ∂h/∂x) + ∂/∂y (K h ∂h/∂y) = S_ya ∂h/∂t

The equation is known as the Boussinesq equation. Note that, in spite of the utilization Dupuit approximation to derive this equation, it is still non-linear since the upper boundary is an unknown priori. Linearization of the equation is possible if the spatial variation of water table is small. Replacing the variable saturated flow with an average thickness b, the linear zed Boussinesq equation become,

∂²h / ∂x² + ∂²h / ∂y² = (S_ya / bK) ∂h/∂t

this is precisely the same form of confined aquifer equation with h representing the height of water table instead of the value of piezometric surface. If accretion W (rain or recharge) is present, this is measured as L/t, the equation become,

∂²h / ∂x² + ∂²h / ∂y² + W/K = (S_ya / bK) ∂h/∂t

Flow in Unconfined aquifers With accretion:

Let the situation be as shown below describing flow between two ditches with recharge source W measured as L/T. The equation describing this flow is,

∂/∂h (Kh ∂h/∂x) + W = 0

Integrating,

∂²h²/∂x² = -2W / K

∂²h²/∂x² = -W/K (x²) +C₁

h² (x) = - W/K (x²) +C₁x + C₂

Applying Boundary conditions,

h = H₁ @ x = 0 and h = H₂ @ x = L

The equation describing steady state in unconfined aquifer with accretion is,

h² (x) = - W/K (x²) + [(H₂² –H₁²)/L + (W/K) L] x + H₁²

To locate the groundwater divide,

dh / dx│_{x = x0} = 0 = -2(W/K) x₀ + (H₁² – H₂²) / L + (W/K) L

x₀ = L/2 – K/2W [(H₁² – H₂²) / L]

To compute The discharge

2h dh/dx = -(2W / K) x + C

Applying the boundary conditions to evaluate C,

C = (H₂² – H₁²) / L – (W/K) L²

Since h dh/dx = -Q/K, then,

Q = Wx + (H₂² – H₁²) / L – (W/K) L²

dq = w dx

Q_R = w ʃ_{X = X0→X = X} dx = x – x₀

= W {x – L/2 +K/2W (H₁² – H₂²)/L}

Q_R = K / 2L (H₁² – H₂²) – W (L/2 – x)

Solution of Groundwater Equation (Radial Symmetry):

It was stated before that the release of water from storage in confined aquifers is due to the compressibility of the aquifer and the compressibility of the water. The drawdown caused by lowering piezometric surface is the vertical distance between the initial and final positions of the piezometric surface at time t. When pumping occurs in an aquifer, non-steady conditions prevail until such a point where the cone of depression reaches a continuous source of replenishment. Theoretically, in the absence of replenishment source, the cone of depression grows with no limit. However, for all practical purposes, the cone will reach a point where no significant drawdown occurs. At this point, the distance to the pumping well is the radius of influence of the aquifer. Thus radius of influence is can be defined as the distance from the pumping well to a point of “zero” drawdown.

When the flow in a confined aquifer is assumed horizontal, the equipotential lines are vertical. The drawdown is a function of x, y and t and not z. The observation well, used to measure the elevation of the piezometric surface, can be located at any depth below the water table. This is not the case in phreatic aquifers, where the observed piezometric head of a point reflects only the head at that point and nowhere else. In both cases, the assumption of horizontal flow near the pumping well cannot be applied because the vertical gradients at this region are high.

Equipotential surface for confined aquifer

(assuming horizontal flow)

Equipotential surface

(Unconfined aquifer)

Convergence of Flow into Pumping Wells:

When pumping begins around unconfined aquifers, water experiences the largest drawdown around the immediate vicinity of the pumping well. The low pressure caused by the pumps creates head deferential between water near the well and the water-bearing formation away from the well causing water to move toward the well. This differential head is the difference between the water level inside the well and the water level at any point in the aquifer. Pumping from confined aquifers does not generally effect the saturation of the aquifer but reduces the hydrostatic head in the aquifer. If desaturation (often referred to mining) occurs due to excess pumping, the aquifer will act as combination of unconfined (around the well) and confined (away from the well). The cone of depression bounded by the hydrostatic head in confined aquifers grows more rapidly with pumping compared to the case in unconfined aquifer. This phenomenon occurs because drainage from unconfined aquifers relates to apparent specific yield, which are many orders of magnitudes larger than the compressibility of water β and the compressibility of soil matrix α, the two terms that define drainage from confined aquifers.

Steady Flow to A Well in Unconfined Aquifer (Radial Symmetry);

Figure below shows a well fully penetrating the saturated unconfined aquifer. The assumption inherent here is that the aquifer is radially symmetric, homogeneous and isotropic.

Cylindrical Cut

Consider flow of groundwater into a ring shaped (cylindrical) element of unconfined aquifer of radius r and height h = H-s, again, the net volume of outflow rate ∂V_w/∂t is,

∂Q_r/∂r ∆r = net volume of outflow rate (=∂V_w/∂t)

where Q = -2π rK (H-s) ∂s/∂r

Substituting for Q, and noting that ∂V_w∂t is -2πr∆rS_y∂s/∂t then,

-2πrS_y∂s/∂t =∂/∂r[-2π rK (H-s) ∂s/∂r]

Or, S_yr∂s/∂t = K∂/∂r[r (H-s) ∂s/∂r]

For linearity, and since s<<H,

r (S_y∂s/∂t) = K H ∂/∂r (r∂s/∂r)

If α = KH/S_y , then r∂s/∂t = α ∂/∂r (r∂s/∂r)

= α r ∂²s/∂r² + α∂s / ∂r. The equation become,

∂s/∂t = α [∂²s/∂r² + (1/r)∂s / ∂r]

This is the linearized form of differential equation for radial flow in unconfined aquifers. If the flow is steady, s is not function of time and the equation become,

∂²s/∂r² + (1/r)∂s / ∂r = 0

This is the radial form of Laplace equation.

Steady Flow State Well Hydraulics

Considering cylinder around the well extending from r_w to R where r_w is the radius of the well and R is the radius of influence (which is the distance from the well to the point where drawdown is negligible), the boundary condition can now be stated as,

h = h_w at r = r_w

h = H at r = R where H is the static head before pumping.

Employing Dupuit Assumption, the flow is,

Q_w = 2πr hK dh/dr

= 2πr K d/dr (h²/2)

Employing the boundary condition h = h_w @ r = r_w and h = H at r = R

H² – h²_w = Q/πK ln (R/r_w)

For small s_w , H² – h²_{w =} (H – h_w) (H + h_w) = s_w (2H) The Solution become,

s_w = Q/2πT ln (R/r_w)

where T is the coefficient of transmissibility KH discussed earlier. The equation only applies to aquifers that are infinite in their lateral extension. To find the drawdown between two points in the aquifer, the equation become,

s₁ – s₂ = Q/2πT ln (r₁/r₂)

The drawdown at any point in the unconfined aquifer is,

s = s_w -Q/2πT ln (r/r_w)

It is apparent from the above equation that, to obtain maximum discharge, the well would have to penetrate to the impermeable boundary. However, this is not economically feasible. Studies found that wells, which penetrate a depth greater than two-thirds of the saturated zone, would not yield significantly more discharge (McWhorter77)

Steady Flow to a Well in Confined Aquifers

Considering a cylinder around the well extending from R_w to R and bounded by the upper and lower impermeable boundaries,

Q_w = 2πr b K dh/dr

After integration, the equation for piezometric head become, b H

H – h = Q/2πbK ln (R/r)

Or, Impermeable

s = Q/2πT ln (R/r)

Relationship of Drawdown to Yield:

Well yield is the amount of water that can be extracted from a given well without producing undesirable effects such as dewatering (mining) of the aquifer. In theory, a 100% efficient well gives 100% yield. However, yield varies depending on topography of the aquifer, the construction of the well bore, and the efficiency of the pump and its operating speed.

From the equation that describes the radial flow in confined aquifers, it is apparent that well yield is directly proportional to the thickness b of the aquifer. However, when pumping becomes excessive, aquifer dewatering will occur. The proportionality of b to yield as prescribed by the equation will no longer hold true due to the reduction of the magnitude of b. On the other hand, pumping from unconfined aquifers dewaters the region around the well and the saturated zone experiences continuous reduction. This, in effect, reduces the amount of water available for pumping and the yield will be affected.

0

100

80

Percent Drawdown

The chart above shows typical relationship between drawdown and yield for a given aquifer. A 50 percent drawdown means lowering water level halfway between the initial saturated zone and the bottom of the aquifer. A 100 percent means lowering the saturated zone to the bottom. If the maximum drawdown of the aquifer is say 10 meters, lowering the saturated zone to 8 meters, (i.e. 20% of the total possible drawdown), will result a yield of 37% of the maximum possible yield. Remember, it is not economically feasible to operate the well at high values of drawdown because of the resulting diminishing return. As the chart indicates, pumping the well at drawdown greater than say 80% of drawdown would require additional pumping power to obtain the remaining 5% yield.

Sustained Yield:

The practical sustained yield is the amount of water that can be extracted from the well without producing irreparable damage to the aquifer. It is generally equal to the amount of available recharge and must not exceed the mean annual recharge. In many arid regions, the rate of replenishment slows down considerably during dry spells. The groundwater resource will eventually be exhausted unless corrective measures are implemented. These include artificial recharge, prevention of waste, reclaiming of used water and importation of water from other sources. Exceeding sustained yield in coastal aquifers is particularly difficult problem. Under natural undisturbed conditions, an equilibrium state between salt water and fresh water exist in the coastal aquifers. If pumping from coastal aquifer exceeds replenishment, The negative gradients produced by excessive pumping lowers the available head at the fresh water zone and an adverse gradient is produced between salt-water zone and the groundwater and the interface will eventually advance inward and contaminates the pumping well. This phenomenon is known as the seawater intrusion.

EXAMPLE

The hydraulic conductivity of an aquifer with r_W = 0.1 meters is 5 meters per day from r_W = 0.1 m to r = 10 m and 15 m / day from r = 10 to r = R = 500 where R is the radius of influence. Find the equivalent hydraulic conductivity of the aquifer and the discharge at the well if the drawdown at the well is 2 and the thickness of the aquifer is b = 20 meters.

In General, h = (Q/2π bK) ln r + C

h₁ – h_W = (Q/2π bK₁) ln (r₁/r_W)

h_R – h₁ = (Q/2π bK₂) ln (R/r₁) Then, h_R – h_W = (Q/2π bЌ) ln (R/r_W)

where Ќ is the equivalent hydraulic conductivity of the whole aquifer. Then,

h_R – h_W = (Q/2π bЌ) ln (R/r_W) = (Q/2π bK₁) ln (r₁/r_W) + (Q/2π bK₂) ln (R/r₁)

(Q/2πbЌ) ln (R/r_W) = (Q/2πb) [(1/K₁ ln (r₁/r_W) + 1/K₂ ln (R/r₁)]

Then, Ќ (the equivalent hydraulic conductivity of the aquifer) is

Ќ = [ ln (R/r_W)] / [(1/K₁ ln (r₁/r_W) + 1/K₂ ln (R/r₁)]

Ќ = [ln (500/0.1] / [(1/5 ln (10/0.1) + 1/15 ln (500/10) = [ln 5000] / (0.2 ln 100 + 0.067 ln 50) = 8.52(0.92 + 0.26) = 10.1 meters per day

The discharge at the well is Q = [2πb Ќ / ln (R/r_w)] (h_R - h_W)

= [(2)(π)(20m)(10.1m/d) / ln (5000)] (2m) = (1269.2/8.52)(2) = 297.93 m³/d

Special Cases:

Case 1: Drawdown in Confined and Unconfined Aquifers in Well Field:

When pumping from well field, the concept superposition can be employed. The concept of superposition is valid only for linear, homogeneous partial differential equations. For a well system comprising of n number of wells, the total drawdown at any point due to the i^th well is,

s_i = Q_i / 2πT ln R/r_i

for r_i< R here, s_i is the drawdown due to well i at any point and R is the radius of influence. The total drawdown is,

s_T = Σ_{i=1 to n} Q_i /2πT ln R/r_i

If the number of wells is large at a relatively small area, the well field can be transformed into uniformly distributed discharge well field of W withdrawal rate (measured in discharge per unit area). The uniform withdrawal rate is,

W = Σ _{i=1 to n} Q_i / A_T

Where A_T is the area containing the system of wells. If we assume that the well field has a circular area of radius r₀ then the differential drawdown at the center of the well system (at r = 0) produced by rate of withdrawal W from a annular ring of area 2πrdr is,

ds₀ = (2π rdr W) / 2πT ln R/r

The drawdown at the center of well field become,

s₀ = W/T∫r ln (R/r) dr integrated from 0 to r₀

s₀ = W r₀² / 2T [ln R/r₀ + ½]

EXAMPLE: (McWhorter 77)

Use the above equation that describe drawdown in well field to estimate the radius of an equivalent well that will produce the same drawdown

The equation describing drawdown for well field is,

s₀ = W r₀² / 2T [ln R/r₀ + ½]

The drawdown in a single well of radius r_W is,

s_W = Q/2πT ln R/r_w

Equating the above two equations,

ln R/r_w = ln R/r₀ + ½. Then r_W = 0.6 r₀

This means that the drawdown ina well field can be replaced by a single well with radius 0.61 of the radius of the well field.

Case 2: Pumping Near Hydro-Geologic Boundaries (Recharge Source)

Consider pumping around fully penetrating stream (recharge boundary). When pumping occurs, the drawdown will increase and the cone of depression will expand until it hits the recharge source. At this point on, the cone of depression cannot spread beyond the recharge source and no drawdown will take place beyond this point. The drawdown at any point in the system is calculated by placing an imaginary recharge (production) well pumping from an infinite aquifer and placed at the exact opposite distance to the recharge source. The recharge well operates simultaneously and at the same pumping rate as the real well. Remember that after this point in time, the flow into the well is no longer radially symmetric because the source of water is the stream.

Impermeable

Recharge source

Cone of

Depression

Due to pumping with

Infinite aquifer (no

Source)

The solution of problems involving pumping near sources in an aquifer can be illustrated as follows:

It is required to find the drawdown at any point in an aquifer for steady flow to a well located at point P(x₀,0). The line x = 0 is a contentious stream boundary infinite in extent at y axis.

The drawdown at any point in the real region in the aquifer P(x,y) is the sum of drawdown of two wells, each operating in a fictitious infinite field. The equation of the drawdown is the sum of the drawdown due to pumping well (+Q) and recharge well (-Q). If r is the distance to the real well and r_i is the distance to the imaginary well, then the drawdown at any point is,

s(x,y) = (Q/2 π T) ln (R/r) + (-Q/2 π T) ln (R/r_i)

= (Q/2 π T) ln (r_i/r)

= (Q/4 π T) ln {[ (x + x₀)² + y²] / [ (x - x₀)² + y²]}

where R is the radius of influence of the well. Note that the assumption here is that R > x₀ otherwise the recharge boundary will have no effect on drawdown. The superposition in side view is illustrated by the first graph in the previous page.

The procedure illustrated in the above example is known as the method of images. Note that the gradient of the head (or drawdown) is zero at any point in the recharge boundary (line x = 0). If more than one recharge boundary exists in the system, the method of superposition still applies and the total drawdown will be the sum of all drawdown due to imaginary wells with distances r_i to the point in question and the drawdown to the pumping well.

s(x,y) = Σ_{i = 1,n} s_i + s_PW

Case 3: Pumping Near Hydro-Geologic Boundaries (Impermeable Boundary)

Consider pumping around fully penetrating impermeable boundary. When pumping occurs, the drawdown will increase and the cone of depression will expand until it hits the impermeable boundary. At this point on, the cone of depression cannot spread beyond the impervious boundary and the rate of drawdown accelerates. The drawdown is calculated by placing a fictitious pumping (discharge) well placed at the exact opposite distance to the impervious boundary. The imaginary discharge well operates simultaneously and at the same pumping rate as the real well. Remember that, after this point in time, and just like the case in recharge boundary, the flow into the well is not radially symmetric flow.

Image well

The solution of problems involving pumping near impermeable boundary in an aquifer can be illustrated as follows:

It is required to find the drawdown at any point in an aquifer for steady flow to a well located at point P(x₀,0). The line x = 0 is a contentious impermeable boundary infinite in extent at y axis.

x

The drawdown at any point in the real region in the aquifer P(x,y) is the sum of drawdown of two wells, each operating in a fictitious infinite field. The equation of the drawdown is the sum of the drawdown due to pumping well (+Q) and recharge well (-Q). If r is the distance to the real well and r_i is the distance to the imaginary well, then the drawdown at any point is,

s(x,y) = (Q/2 π T) ln (R/r) + (Q/2 π T) ln (R/r_i) = (Q/2 π T) ln (R²/rr_i)

= (Q/2 π T) ln { R²/ {[ (x - x₀)² + y²]^1/2 [ (x + x₀)² + y²]^1/2}

Where R is the radius of influence of the well. Note that the assumption here is that R > x₀ otherwise the impermeable boundary will have no effect on drawdown. The superposition in side view is illustrated by the first graph in the previous page. Writing the above equation in terms of h_r and h_w,

h_r – h_w = (Q/2 π T) ln (r/r_w) + (Q/2 π T) ln (r_i/r_w)

= (Q/2 π T) ln (rr_i/r_w²)

And for r ≈ r_i

h_r – h_w = (Q/ π T) ln (r²/r_w²)^1/2 then

h_r at any point = h_w + (Q/ π T) ln (r/r_w) which is twice the drawdown produced by single well from infinite aquifer.

Case 4: Image Well System:

Method of images also applicable in a system of wells that is bounded by two types of boundaries at angles less than 180 degrees. Some are shown in the following illustrations,

Drawdown at any point in the real domain is,

S_P = s_r + s₁ – s₂ – s₃

= (Q/2πT) [ln ( R/r) + ln (R/r_i1) – ln (R/r_i2) –ln (R/r_i3)

= (Q/2πT) [ln ( r_i2r_i3 / rr_i1)

= (Q/2πT) {ln [(b+x)² + (y+a)²]^1/2 [(b+x)² + (y-a)²]^1/2 / [(x-b)² + (a-y)²]^1/2 + [(x-b)² + (y+a)²]^1/2

The equation becomes,

S_P (x,y) = (Q/4πT) ln { [(b+x)² + (y+a)²] [(b+x)² + (y-a)²] / [(x-b)² + (a-y)²] + [(x-b)² + (y+a)²]}

Case 5: Steady State Drawdown from Leaky Confined Aquifers

The equation describing distribution of drawdown (or head) in leaky, confined aquifers is,

1/r ∂/∂r (r∂h/∂r) + (H – h) / λ² = 0, where,

λ² = b b^’ K / K^’

λ is called the leakage factor and b^’ is the thickness of the leaking formation and H is the head at radius of influence. In terms of s (r) the equation is,

∂²s/∂r² + (1/r) ∂s/∂r + s / λ² = 0

The solution of the equation is,

s_r = Q/2πT [K₀ (r/λ) ]

Conf.

Aquifer

b

where K₀ is the modified Bessel function of the second kind order zero. The inherent assumption for the above solution is that aquifer is infinite in extent and r_w/λ << 1. For the vicinity of the well, the drawdown becomes,

s_r = Q/2πT ln (1.123λ/r)

Modified Bessel Function K₀ (r/λ):

N:	N x 10^-3	N x 10^-2	N x 10^-1	N
1.0	7.0237	4.7212	2.4271	0.4210
1.5	6.6182	4.3159	2.0300	0.2138
2.0	6.3305	4.0285	1.7527	0.1139
2.5	6.1074	3.8056	1.5415	0.0623
3.0	5.9251	3.6235	1.3725	0.0347
3.5	5.7709	3.4697	1.2327	0.0196
4.0	5.6374	3.3365	1.1145	0.0112
4.5	5.5196	3.2192	1.0129	0.0064
5.0	5.4143	3.1142	0.9244	0.0037
5.5	5.3190	3.0195	0.8466
6.0	5.2320	2.9329	0.7775	0.0012
6.5	5.152	2.8534	0.7159
7.0	5.0779	2.7798	0.6605	0.0004
7.5	5.0089	2.7114	0.6106
8.0	4.9443	2.6475	0.5653

Transient Well Hydraulics

Groundwater flow in confined and unconfined aquifers is transient (variable with time) when the piezometric surface or water table position changes with time. The solution is obtained by solving the linearized form of Boussinesq equation, which employs Dupuit assumption. In the following analysis, the aquifer is assumed homogeneous, isotropic, for the case of confined aquifers, the thickness is constant, and in both cases, storativity of the aquifer is constant. Furthermore, the release of water from the aquifer is immediate upon decline of head. It is important to note that, upon employment of Boussinesq equation, it is apparent that the restrictive nature of Boussinesq equation is more pronounced in unconfined aquifers than in confined aquifers.

The linearized form of differential equation of flow of groundwater in radial symmetry was presented earlier (page 28) as,

∂s/∂t = α [∂²s/∂r² + (1/r)∂s / ∂r]

Where

α = T/S for confined aquifers and T/S_y for unconfined aquifers.

Employing boundary conditions: s = 0 at t = 0, r > 0

r∂s/∂r_r→0 = -Q/2πKb at r → 0 for confined aquifers, and,

r∂s/∂r_r→00 = -Q/2πKH for unconfined aquifers (remember that, for small s, we assumed that H ≈ H – s.

Using the “Boltzman” transformation u

u = r² / 4αt

the equation become,

d²s/du² + (1 + 1/u) ds /du = 0

The boundary conditions are transformed in term of u as: s(t=0) or s(∞) = 0

And: u ds/du_u→0 = -Q/2πT

Integrating yield,

u ds / du =C₁ e^-u where C₁ is a constant of integration. Applying the second boundary condition,

ds / du = - (Q/2πT) e^-u / u

Integrating and employing the first condition, the solution becomes,

s = Q/4πT ʃ _U→∞ e^-X / x dx

Where x is a dummy variable of integration.

The integral above is known as the exponential integral and in groundwater literature it is known as the Theis well function W (u) (after Theis 1935)

The equation of drawdown in transient flow conditions becomes,

s = Q/4πT ʃ_U→∞ e^-U / u du and u = r² S/4Tt

Or, in short s = Q/4πT W (u)

which is the non-equilibrium equation describing transient flow of groundwater to a well developed by Theis (1935).

The well function W (u) is can be expanded in an infinite series as follows,

W (u) = -0.5772 – ln u + u – u²/2.2! +u³/3.3! – u⁴/4.4! + ……………………

For small values of u (say u<0.01, i.e. for small r and/or large time), the W function can be approximated by the first two terms and the equation can be approximated as,

s (r,t) = Q/4πT [ -0.5772 + ln(1/u)]

Substituting for u = r² S/4Tt, the approximate equation become,

s (r,t) = Q/4πT ln (2.246Tt / r²S)

The above equation approximates the drawdown in both confined and unconfined aquifers with S represents storage coefficient for confined aquifers and apparent specific yield for unconfined aquifers. In addition, T would be = bK for confined aquifers and = ĤK for unconfined aquifers where Ĥ represents an “average” saturated thickness.

Observations:

• The well function equation developed above can be applied to unconfined aquifers in the most restrictive sense. The assumptions inherent in the application of this equation to unconfined aquifers are that gravity drainage is instantaneous and no delayed yield occurs. The other assumption is that flow is horizontal and drawdown is too small compared to the overall saturated thickness of the unconfined aquifer. In this approach, which is by far the simplest, is to use the same equation as the confined aquifer equation but with different arguments. For instance, the transmissibility Kb is replaced by the value KH where H either is the initial saturated thickness or “averaged” thickness of unconfined saturated zone.

• The exact solution of the integral equation predicts that the cone of depression around the well develops instantaneously and extends infinitely.

• Since u is a function of r, then as u → ∞, r → ∞ and W (u) → 0. This makes r → 0 and for all practical purposes, the drawdown becomes negligible for a finite radius. Mathematically, however, the equation suggests that the radius of influence R grows infinitely with time.

• Since R is function of t^1/2, the cone of influence expands rapidly initially and slows down thereafter.

• The equation of R tells us that, given the same pumping conditions the cone of influence is larger in aquifers with small S compared to those of large S. This means that in an unconfined aquifer (which has S_y>>S), the cone of influence does not have to expand as fast as the cone of influence in confined aquifers to account for the same value of Q.

• Physical limit of cone of influence is replenishment source (such as rivers, lakes or sea).

• Mathematically, steady state conditions will never prevail because u is not constant.

• A condition known as pseudo-steady state occur by setting u equal or less than an arbitrary value of 0.01 if attention is restricted to regions around the well and pumping takes place at sufficiently long time. Remember that pseudo-steady state does not imply steady state conditions where s does not vary with time, it rather imply that, at regions not far from the well, the rate of fall of the piezometric surface or the water table is the same everywhere at this region and the rate of change in drawdown is independent of distance r.

The Transient Groundwater Equation:

s = Q/4πT [W (u)]

W (u) = -0.5772 – ln u + u – u²/2.2! +u³/3.3! – u⁴/4.4! + ……………………

u = r² S/4Tt

Values of W (u) (After Wenzel 1942):

1.00E-15	33.96	9.00E-12	24.86	1.00E-07	15.24	9.00E-04	6.44
2.00E-15	33.27	1.00E-11	24.75	2.00E-07	14.85	1.00E-03	6.33
3.00E-15	32.86	2.00E-11	24.06	3.00E-07	14.44	2.00E-03	5.64
4.00E-15	32.58	4.00E-11	23.36	4.00E-07	14.15	3.00E-03	5.23
5.00E-15	32.35	5.00E-11	23.14	5.00E-07	13.93	4.00E-03	4.95
6.00E-15	32.17	6.00E-11	22.96	6.00E-07	13.75	5.00E-03	4.73
7.00E-15	32.02	7.00E-11	22.81	7.00E-07	13.6	6.00E-03	4.54
8.00E-15	31.88	8.00E-11	22.67	8.00E-07	13.46	7.00E-03	4.39
9.00E-15	31.76	9.00E-11	22.55	9.00E-07	13.34	8.00E-03	4.26
1.00E-14	31.66	1.00E-10	22.45	1.00E-06	13.24	9.00E-03	4.14
2.00E-14	30.97	3.00E-10	21.35	2.00E-06	12.55	1.00E-02	4.04
3.00E-14	30.56	4.00E-10	21.06	3.00E-06	12.14	2.00E-02	3.35
4.00E-14	30.27	5.00E-10	20.84	4.00E-06	11.85	3.00E-02	2.96
5.00E-14	30.05	6.00E-10	20.66	5.00E-06	11.63	4.00E-02	2.68
6.00E-14	29.87	7.00E-10	20.5	6.00E-06	11.45	5.00E-02	2.47
7.00E-14	29.71	8.00E-10	20.37	7.00E-06	11.29	6.00E-02	2.3
8.00E-14	29.58	9.00E-10	20.25	8.00E-06	11.16	7.00E-02	2.15
9.00E-14	29.46	1.00E-09	20.15	9.00E-06	11.04	8.00E-02	2.03
1.00E-13	29.36	2.00E-09	19.45	1.00E-05	10.94	9.00E-02	1.92
2.00E-13	28.66	3.00E-09	19.05	2.00E-05	10.24	1.00E-01	1.82
3.00E-13	28.26	4.00E-09	18.76	3.00E-05	9.84	2.00E-01	1.22
4.00E-13	27.97	5.00E-09	18.54	4.00E-05	9.55	3.00E-01	0.91
5.00E-13	27.75	6.00E-09	18.35	5.00E-05	9.33	4.00E-01	0.7
6.00E-13	27.56	7.00E-09	18.2	6.00E-05	9.14	5.00E-01	0.56
7.00E-13	27.41	8.00E-09	18.07	7.00E-05	8.99	6.00E-01	0.45
8.00E-13	27.28	9.00E-09	17.95	8.00E-05	8.86	7.00E-01	0.37
9.00E-13	27.16	1.00E-08	17.84	9.00E-05	8.74	8.00E-01	0.31
1.00E-12	27.05	2.00E-08	17.15	1.00E-04	8.63	9.00E-01	0.26
2.00E-12	26.36	3.00E-08	16.74	2.00E-04	7.94	1.00E+00	0.219
3.00E-12	25.96	4.00E-08	16.46	3.00E-04	7.53	2.00E+00	0.049
4.00E-12	25.67	5.00E-08	16.23	4.00E-04	7.25	3.00E+00	0.013
5.00E-12	25.44	6.00E-08	16.05	5.00E-04	7.02	4.00E+00	0.0038
6.00E-12	25.26	7.00E-08	15.9	6.00E-04	6.84	5.00E+00	0.0011
7.00E-12	25.11	8.00E-08	15.76	7.00E-04	6.69	6.00E+00	0.00036
8.00E-12	24.97	9.00E-08	15.65	8.00E-04	6.55	7.00E+00	0.00012

EXAMPLE:

Calculate the limit of the pseudo-steady state region around the well if Q = 100 m³ / hr, S_ya = 0.09, T = 36 m² / hr and t = 10 hrs.

r² / 4αt = r² / 4 (T/S_ya)t

= 0.01 then

r = [(0.01)(4)(36/0.09)(10)]^1/2

= 12.65 m.

The Pseudo-steady state does not apply beyond this limit.

EXAMPLE

Compare the effective radii of influence at t=10 hrs for wells pumping at 100 m³/hr in aquifers for which a) T = 36 m²/hr, S_ya = 0.1, b) T = 360 m²/hr, S_ya = 0.1, c) T = 36 m²/h, S = 0.01. Let the radius of influence be the radius at which s_e = 0.01

W(u) = 4πTs_e / Q = 4π(36)(0.01) / 100 = 0.045 from table @ w(u) = 0.045, u =2.5684 then,

r_e = [ u(4)Tt / S_ya]^1/2 = [2.5684 (4)(36)(10)/0.1]^1/2 = 192.3 m

W(u) = 4πTs_e / Q = 4π(360)(0.01) / 100 = 0.45 from table @ w(u) = 45, u = 0.6253 then,

r_e = [ u(4)Tt / S_ya]^1/2 = [0.6253 (4)(360)(10)/0.1]^1/2 = 300.1 m

W(u) = 4πTs_e / Q = 4π(36)(0.01) / 100 = 0.045 from table @ w(u) = 0.045, u = 2.5684then,

r_e = [ u(4)Tt / S_ya]^1/2 = [2.5684 (4)(36)(10)/0.01]^1/2 = 608.2 m

From the above analysis, it is evident that the radius of influence is more sensitive to the storage characteristics of the aquifer than the transmissibility.

EXAMPLE:

After 2 hours of pumping in an aquifer, it was observed that the radius at which drawdown is negligible is 400 meters. At what radius would the drawdown be negligible after 5 hours of pumping? Assume pseudo-steady state prevail.

At pseudo-steady state,

s = Q/4πT [ ln (2.246 Tt / R²S)

0 = ln [(2.246 t / R²) (T/S)]

1 = (2.246) (2)/160000) T/S

T/S = 2.81 x 10^-5

For 5 hrs of pumping,

0 = ln [2.246 (5) / R²] (2.81 x 10^-5)

Then,

R² = (5) (2.246) / (2.81 x 10^-5)

R = 632.46 m

EXAMPLE:

A pumping well is bounded by two straight parallel streams 20 meters apart. The well is located at the middle. Pumping starts at Q = 1000 m³/d. If T = 500 m³/d , find the drawdown at the well after 10 days of pumping. Assume R = 500 m and r_PW = 0.5 m and S_ya = 0.1.

RW_i

Impermeable

The transient drawdown equation is,

s = Q/4πT W(u), where,

u = r²S / 4Tt,

Taking four image wells from each side,

s_T = s_PW + (- s_iRW20 - s_iRW40 + s_iPW60 + s_iPW80 )_LEFT + ( s_iPW20 - _iRW40 – _iRW60 + s_iPW80 )_RIGHT

= s_PW - 2 s_iRW40 + 2 s_iPW80, the rest cancel each other.

The drawdown at the well is,

s_T = Q/4πT [ W(u_PW) - 2 W(u_iRW40) + 2 W(u_iRW80)]

u_PW = (0.5)² (0.1) / (4)(500)(10) =0.00000125,

u_iPW40 =(40)² (0.1 / (4) (500)(10) = 0.0008

u_iPW80 =(80)² (0.1) / (4) (500)(10) = 0.032

The total drawdown become,

s_T = (1000) / (4 π) (500) [13.03555 – (2) (6.5545) + (2) (2.8845)]

= 0.1592 [13.03555-13.109+5.769]

Drawdown at the well after 100 days of pumping:

s_T = 0.91 m

EXAMPLE:

Rework the same problem above with only one stream located 50 meters from a pumping well. The discharge is Q = 1000 m³/d, T = 500 m²/d, , r_PW = 0.2 m and S_ya = 0.1. Find the drawdown after 0.1, 1.0, 10, 100 and 1000 hrs of pumping. Analyze the results.

The transient drawdown equation is:

s = Q/4πT W (u), where, u = r²S / 4Tt,

The total drawdown is: s_T = s_PW0.1 + s_PW1 + s_{PW10 +} s_PW100 – s_iRW0.1 - s_iRW1 - s_iRW10 - s_iRW100

u_PW0.1 = (0.2)²(0.1)/4(500)(0.1) =0.000002, u_iRW0.1 =(100)²(0.1)/4(500)(0.1) =5

u_PW1 = (0.2)²(0.1)/4(500)(1) =0.000002, u_iRW1 =(100)²(0.1)/4(500)(1) =0.5

u_PW10 = (0.2)²(0.1)/4(500)(10) =0.0000002, u_iRW10 =(100)²(0.1)/4(500)(10) =0.05

u_PW100 = (0.2)²(0.1)/4(500)(100) =0.00000002, u_iRW10 =(100)²(0.1)/4(500)(100) =0.005

Drawdown at the well after 0.1 day of pumping:

s_T0.1 = 1000 / (4π)(500) [W(0.00002) – W(5] = 0.1592 [10.2426 – 0.001148] = 1.630 m

Drawdown at the well after 1.0 day of pumping:

s_T1.0 = 0.159 [W (0.000002) –W (0.5) ] = 0.1592 [12.5451 –(0.5598)] = 1.908 m

Drawdown at the well after 10 days of pumping:

s_T10 = 0.159 [W (0.0000002) –W (0.05) ] = 0.1592 [14.8477 –(2.4679)] = 1.971

Drawdown at the well after 100 days of pumping:

s_T100 = 0.159 [W (0.00000002) – (2) W (0.005) ] = 0.1592 [17.1503 – (4.7261)] = 1.9778 m

Steady state drawdown:

S_ss = Q / 2πT ln (r_i / r_w) = 1000 /(2π)(500) ln (100/.2) = 1.9781 m

Which indicates that at large time, the infiltration from the stream balances the discharge from the well and the drawdown remain constant.

Analysis:

From the example above, the drawdown at the production well near a stream is calculated using the non-equilibrium equation of the production well and its counterpart placed in the opposite side of the stream in an infinite aquifer system,

s = Q/4πT [W (u_PW) – W (u_iRW)]

The first term of the equation represents the drawdown that is due to the production well in the “real” field and the second term represents the “recovery” of drawdown due to the recharge well. Initially, all of the pumped water is derived from the aquifer storage until the time when the cone of influence intersects the recharge source. When the cone reaches the recharge source, a hydraulic gradient develops between stream and groundwater because of drawdown at the well. The build-up of drawdown induces seepage from the stream at increasingly larger rate. The increase of water due to replenishment affects the drawdown at the well causing the rate of drawdown at the well to decrease continuously as was shown in the example. The rate of decrease of drawdown continues until steady state conditions prevail and the infiltration from the stream balances the discharge at the well.

Plotting the result of the example on semi logarithmic paper shows that the cone of influence was enlarging during the first day. The second part of the curve demonstrates the influence of the recharge source at subsequent days, which is evident in the flattened part of the curve.

Time (day)-log

0.01 0.1 1.0 10 100 1000

Drawdown as a function of time in presence of recharge source

Delayed Yield:

This approach was pioneered by Boulton 1954 and advanced by Neuman 1972. The concept of delayed yield applies in unconfined aquifers and will be discussed subsequently.

Special Cases:

Case 1: Transient Flow in Unconfined Aquifers Accounting for Delayed Yield:

When water is pumped from unconfined aquifer, a cone of depression develops in the saturated water zone. This zone delivers water to the well through three subsequent and distinct mechanisms. The first is an immediate response to pumping which resembles the mechanism of delivery in confined aquifers, namely the compaction of aquifer and expansion of water (segment A in the curve). The value of the storage coefficient computed from this segment is well within the confined range and is not representative of the storage coefficient of the aquifer. This value is many orders of magnitudes smaller than the apparent specific yield of the aquifer. This “elastic” response lasts few moments following the start of pumping. In the second segment, gravity drainage takes effect and water delivered to the well is due to dewatering caused by falling water table. Under this segment, which is not unlike that of vertical leakage into unconfined aquifers, the cone of depression slows in its rate of expansion as water released from storage by gravity feeds the cone (segment B). Finally, the rate of expansion of the cone of depression increases as gravity drainage keeps pace with declining water level (segment C). This segment conforms to a Theis-type curve where vertical gradients of piezometric heads are small in comparison to horizontal gradients. Furthermore, the rate of increase of drawdown is somewhat greater than that of the second segment but far smaller than the first segment. The measured coefficient of storage increases at diminishing rate but the existence of some vertical flow would not permit its value to reach the specific yield of the aquifer.

Aquifer elasticity zone

Drawdown

Time

Gravity drainage zone

Following pumping, the amount of water delivered from storage in unconfined aquifer due to an increment of drawdown ∆s between t and t+∆t consists of two components: a) water instantly released from storage due to compaction of aquifer and expansion of water, and b) water released from storage due to the action of delayed yield. The equation describing the flow, accounting for gravity drainage is (Streltsova 1972),

∂²s/∂r² + (1/r) ∂s / ∂r = (S/T∂s/∂t + S_y/T ∂s⁰/∂t)

Where T is average transmissivity, s is the drawdown due to change in piezometric head and s⁰ is the drawdown due to vertical drainage. The vertical Darcy’s velocities for both drawdown mechanisms are,

q^z = - K (s – s⁰) / b_a where b_a is 1/3 the initial saturated thickness, the Darcy’s velocity of the declining water table is,

q_z = -S_y ∂s⁰ /∂t

Combining the three equations yields the approximate equation describing the transient unconfined flow accounting for delayed yield,

∂²s/∂r² + (1/r) ∂s / ∂r = S/T∂s/∂t + (s – s⁰) / λ²

If the delayed yield is not accounted for, the above equation would describe leaky aquifer discussed below.

Case 2: Transient Analysis for Leaky Aquifers:

The assumption inherent in Theis analysis is that the geological formations overlying and/or underlying the aquifer are impermeable. Should theses formations produce significant inflow to the main aquifer, Theis solution would have to be modified to account for the flow from these formations.

Assuming that the well fully penetrates the main aquifer and the aquifer is infinite in areal extent; the amount of discharge at the well comes from two sources: the storage from the aquifer and the leakage from the aquitard. The rate of leakage is proportional to the difference between the piezometric surface and the water table of the source bed. The differential equation describing the drawdown is,

∂²s/∂r² + (1/r)∂s / ∂r – s/λ² = S/T ∂s/∂t

Where λ is the leakage factor which was defined earlier. The boundary conditions are:

s → 0 as r → ∞ for t > or = 0, and lim_r→0 r ∂s∂ / ∂r = -Q/2πT

The initial conditions,

s = 0 for t < 0 for all values of r. The solution of the equation (Hantush and Jacob 1955) is,

s = Q/4πT ò_u→∞ 1/y exp (-y - r² / 4B² y) dy, for which the solution is,

s = Q/4πT W (u, r/λ) u = r²S / 4Tt

where W (u,r/λ) is known as the well function for leaky confined aquifers. Note that when the discharge at the well is balanced by the leakage from aquitard, water level stabilizes and steady state will develop. The equation become, ∂²s/∂r² + (1/r) ∂s / ∂r – s/λ² = 0 which is the steady state equation of discussed earlier.

Values of W(u,r/λ) for Leaky Confined Aquifer (Hantush 1956):

r/λ:

U	0.01	0.05	0.20	0.50	1.00	1.50	2.00	2.50
5x10^-6	9.441
1x10^-5	9.418
5x10^-5	8.8831
1x10^-4	8.398	6.228
5x10^-4	6.975	6.082
1x10^-3	6.307	5.796	3.505
5x10^-3	4.721	4.608	3.457
1x10^-2	4.036	3.980	3.288	1.849
5x10^-2	2.468	2.458	2.311	1.493	0.841
1x10^-1	1.823	1.818	1.753	1.442	0.819	0.427	0.228
5x10^-01	0.560	0.559	0.553	0.521	0.421	0.301	0.194	0.117
1x10^-0	0.219	0.219	0.218	0.210	0.186	0.151	0.114	0.080
5x10^-0	0.001	0.001	0.001	0.001	0.001	0.001	0.001	0.001

Case 3: Drawdown with Variable Pumping Rates:

The linearity principle of groundwater equation allows for the application of the principle of superposition of drawdown resulting from variation of pumping rates. The principle of superposition also applies to drawdown recovery following well shutdowns. If amount of water pumped changes from Q₁ to Q₂, the resulting drawdown is described by the following equation,

s = (Q₁/4πT) W(r²/4Tt) + (Q₂ – Q₁) W [r²S/4T(t-t_i)] for t>t_i

or, in genera_l’

s = 1/4πT Σ ∆Q_ii+1 W [r²S/4T(t-t_i)]

The equation states that the resulting drawdown is equal to the drawdown from t = 0 to t = t₂ with Q = Q₁ plus drawdown resulting from an “imaginary” well pumping at a rate of Q = ∆Q placed in the same position and pumping from t = t₁ to t = t₂.

t₄

DD due Q₁ @ t’ DD due Q₂ – Q₁@ t” DD at t^’’’ Recovery

If Q₂ is zero and Sr² / 4t T < 0.01 then

s = (Q/4πT ln [t /(t-t₁)]

Recovery

Buildup

Case 4: Pumping Near A Stream:

It was stated before that when water is pumped near a stream, part of the water is derived from the aquifr storage and the other part is from the nearby stream. The discharge from the stream can be calculated from integrating Darc’s formula,

Q_s = -Kb ʃ_-∞→∞ ∂s / ∂x∣_x=0 dy

The result is,

Qs = Q [1 – erf (x₀ / √4KH/S)t]

Where x₀ is the perpendicular distance to the stream and the error function erf is,

erf (x) = 2/√π ʃ_0→x exp (-y2) dy

Note that as the argument x becomes small (large time), the value of erf (x) gets small and the discharge from the stream to he aquifer approaches the well discharge Q.

Values of Error Function erf (x):

X erf (X)	0.00	0.20	0.40	0.60	0.80	0.90	1.00	1.10	1.20	1.30
	0.000	0.223	0.428	0.604	0.742	0.797	0.843	0.880	0.910	0.934

1.40	1.50	1.60	1.70	1.80	1.90	2.00	∞
0.952	0.966	0.976	0.984	0.989	0.993	0.995	1.00

EXAMPLE:

A well pumping at a rate of 2000 m³ / d located 100 m from a recharge source. The well shuts down after a week of pumping. If the aquifer transmissibility is 1000 m² / d and the apparent specific yield is 0.1. Compute the discharge portion coming from the stream after five days and nine days from the start of pumping.

X₀ / √4(T/S) t = 100 / [(4) ( 1000 /(0.1)(5)]^1/2 = 0.224

erf (0.224) = 0.2296

Q_s = (2000) (1-0.2296) =1540.8m³/d

Following pump shutdown, the water recovered from the stream is calculated by assuming the well to be continually discharging at the original rate from t = 0 to t = 10 days plus a second recharging (negative discharging) well placed at the same location pumping at the same rate but starting a week later.

Q

Q_s = Q [ 1 – erf X₀ / √4(T/S)t] – Q [ 1 - erf X₀ / √4(T/S)(t-t₂)

-Q

where t = 9 days and t₂ = 7 days

Q_s

Q

Q_s = Q{ [ erf x_0=/√4(T/S)(t-t₂)] - erf X₀ / √4(T/S)t]}

Q_s = 2000{ [ erf [100 /√4(10000)(3) ] – [erf [100 /√4(10000)(9)]}

= 2000 (0.296 – 0.185) = 222 m³ /

Aquifer Tests

Aquifer tests are field tests that are performed to determine field data under controlled conditions. The outcome of the field tests are the hydrologic parameters of the aquifer such as storage coefficient, apparent specific yield, hydraulic conductivity and transmissibility of the aquifer. The obtained values of the tested aquifer can be used for designing well field and predicting future drawdown. They can also be used for assessing groundwater supply, and estimating inflow and outflow to and from groundwater basins.

The outcome of aquifer properties obtained from field tests will eventually be compared to the values obtained from theoretical considerations. For this reason, it is important that elements such as boundary conditions, initial conditions and the physical characteristics of the chosen sites match closely the ones assumed in theory. The other important consideration is to minimize uncertainties associated with the selection of pumping wells and the placement of observation wells. For example, data obtained from partially penetrating wells are difficult and uncertain to analyze than the fully penetrating wells. In addition, the proper choice of the number and the location of observation wells would minimize uncertainties associated with measuring hydraulic parameters of the aquifer.

In order to ensure that the water level properly represents the average piezometric head below water table, the casing in the observation wells should be perforated and should penetrate the entire saturated zone. In confined aquifers, piezometer should be sealed properly in order to ensure against water transport from one stratum to another.

During aquifer test analysis, the number and the appropriate spacing of observation wells play an important part to insure the integrity of data collected. According to Walton (1987), no less than three observation wells should be selected and spaced at one logarithmic cycle from one another. Because the radius of influence expands more rapidly in confined aquifers than in unconfined aquifers, observation wells should be placed at greater distances from each other than in unconfined aquifers.

When hydrogeologic boundary is present, and in order to minimize their effect, pumping well should be placed at least one saturated thickness away from the boundary. In addition, the observation wells should be spaced along a line through the production well and parallel to the boundary. This is made to minimize the effect of the boundary on distance-drawdown data.

Finally, the hydrologic data collected for unconfined aquifers using aquifer test analysis represent mean values of that portion of aquifer bounded by the initial saturated zone and the of cone of depression and does not represent the entire saturated thickness. Furthermore, the best estimates of the aquifer properties can be obtained from tests that are conducted over a long time when delayed drainage is not a factor influencing the apparent specific yield.

Methods of Analyzing Aquifer Test Data:

a) Analysis of Aquifer Test Data Using Theis Solution:

This method utilizes type curve matching technique. This method involves matching logarithmic time-drawdown or distance-drawdown graphs obtained from the field measured values of time, drawdown or distance, drawdown data with the theoretical logarithmic well function curves (type curves) of W(u) vs. u. The procedure is illustrated as follows:

1. At a given time or at a given location obtain measured field data of drawdown and plot drawdown s vs. r² / t on a log-log paper.
2. Plot “type” curve of W(u) versus u choosing appropriate cycle(s) on a log-log paper of the same size.
3. Super-impose the two plots data by keeping the coordinate axis of the two curves parallel until best fit is obtained.
4. Select any arbitrary point in the graph. Read the values of W(u)^*, u^* and the corresponding values of s^* and (r²/t)^*
5. The values of T and S are calculated using the formulas: T = Q W(u)^*/4πs^* and S = 4Tu^*/(r²/t)^*

u

b) Analysis of Aquifer Test Data Using Cooper-Jacobs Method

It was noted by Cooper and Jacob (1946) that for small u (large values of t or small values of r), the non-steady equation for s become,

s (r,t) = Q/4πT (0.5772 – ln r²S/Tt)

or, in terms of decimal log scale, the equation become,

s (r,t) = 2.303 Q/4πT log (2.246Tt / r²S)

Plotting the drawdown s versus the logarithm of time t would form a straight that has the intercept 2.303 Q/4πT

and the slope 2.246Tt / r²S. The procedure is summarized as follows:

1. For stationary reading, measure s at different times.

2. Plot s vs. log t on a semi-log paper.

3. Find the slope by reading the difference between two points across one log cycle ∆s^*.

4. Find the horizontal intercept t₀^* on the log t axis.

5. The values of S and T are calculated from the formulas: S = 2.246Tt₀^* / r², and T = 2.303 Q/4π ∆s^*.

An alternative method is to measure s at different observation wells (different r), plot s versus r on logarithmic paper, find the slope across one cycle (∆s^*) and find the intercept r₀^* on the log r axis. The values of S and T are calculated from the formulas: S = 2.246Tt / r₀² and T = 2.303 Q/4π∆s^*.

Note that s versus t is a straight line as long as the assumption of small value of u is still valid. However, it is apparent that when time is small, this assumption can no longer be valid as shown in the above chart. It is also important to note that using non-equilibrium equation is valid for unconfined aquifers as well as long as the measured drawdown is small in relation to the overall saturated thickness of the unconfined aquifer.

In order to Illustrate the use of the methods discussed above, the following is a data collected from an observation well during a test in unconfined aquifer in Fort Collins, Colorado (McWhorter & Sunada 77):

S (m):	0.025	0.050	0.055	0.110	0.170	0.180	0.220	0.300	0.370	0.450	0.530	0.620	0.640
r^2/t	88.9	53.3	47.1	25.0	16.7	15.1	11.1	6.25	4.12	2.47	1.55	0.98	0.82

Plotting W(u) versus u, the coordinates of match point W(u) = 1.0 u = 0.1 s = 0.183 and r² / t = 6.2

Then T = (1.872)(1.0) / 4 π (0.183) = 0.814 m²/min and S_ya = (4) (0.814) (0.1)/ 6.2 = 0.053

0.1

10.0

1.0

Another data

Data collected during a test of a confined aquifer by U.S. Geological Survey is as follows:

t (min) s (meter)	1.0 0.200	2.0 0.300	3.0 0.370	4.0 0.415	5.0 0.450	6.0 0.485	8.0 0.530	10.0 0.570	12.0 0.600	14.0 0.635	18.0 0.670	24.0 0.720
t (min) s (meter)	30.0 0.760	40.0 0.810	50.0 0.850	60.0 0.875	80.0 0.925	100.0 0.965	120.0 1.000	150.0 1.045	180.0 1.070	210.0 1.100	240.0 1.200

r = 61 m and Q = 1.894 m³ /min

Plotting the given data on a semi log paper,

∆s = 0.4 m

Per Log cycle

The slope of the line is ∆s = 0.4, then

T = (2.303) (1.894) / 4π (0.4) = 0.868 m² / min

The extrapolated t at s = 0 is 0.4 minute, then

S = (2..246) (0.868) (0.4) / (61)² = 2 x 10^-4

c) Determination of Transmissibility from Recovery Tests:

It is possible to determine the aquifer transmissibility from the recovery of water following pump shutdown if u is sufficiently small. This method is particularly valuable when conditions do not allow for constructing observation wells and drawdown are not disturbed by the action of the pump. The equation used for aquifer recovery is (Page 46),

s = (Q/4πT ln [t /(t-t^’)] where t’ is the time the pump shuts down

or, in terms of decimal log scale, the equation become,

s = 2.303 Q / 4πT log [t/(t-t^’]

This equation is a straight line when plotted in semi-logarithmic paper. When the residual drawdown is plotted against log (t/(t-t^’)], the slope ∆s over one log cycle allows for the determination of T,

T = 2.303 Q / 4π∆s

The method is illustrated by the following example (adopted from Todd 1980):

A well Pumping at uniform rate of 2500 m³/day was shut down after 240 minutes; thereafter, measurements of s and (t-t^’) were made in an observation well and tabulated blow. Find the transmissibility of the aquifer.

t (m)	241	242	243	245	247	250	255	260	270	280	300	320	340	380	420
t - t^’	1.0	2.0	3.0	5.0	7.0	10.0	15.0	20.0	30.0	40.0	60.0	80.0	100.0	140.0	180.0
t/(t-t^’)	241	121	81	49	35	25	17	13	9	7	5	4	3.4	2.7	2.3
s (m)	0.89	0.81	0.76	0.68	0.64	0.56	0.49	0.55	0.38	0.34	0.28	0.24	0.21	0.17	0.14

t/(t-t^’)

From the graph, ∆s = 0.4, then’

T = 2.303 Q /4π∆s = 2.303 (2500) / 4π (0.4)= 1145.4 m² / d

d) Use of Slug Tests to determine Transmissibility:

Aquifer tests have many practical disadvantages, some include:

a) Expensive because they require the construction of observation wells and require person-hour to monitor them.

b) In contaminated sites, water removed from the well would pose a disposal problem.

c) Water level at the well is often disturbed by the action of the pump.

d) Pumping from low permeability aquifer is often difficult to conduct.

Slug tests are field tests used to determine the transmissibility in aquifers with medium to low hydraulic conductivities. They are particularly useful in hazardous sites where pumping wells are not present and water could not be disposed. The method consists of instantaneous injection (or withdrawing) of water from the well. If injection method is used, a small amount of water is injected into the well, or as an alternative, a small slug is placed just below the initial water table. This is equivalent to injecting a volume of water equal to the volume of the slug. If u is sufficiently small, the transmissibility of the aquifer is calculated from measuring the dissipation with time of the mound created by the slug. If the volume of the slug is V, the difference between the elevations before and after the mound build up is defined by the following equation (Ferris and Knowles 1963),

s = - V / 4πTt exp [-r²S / 4Tt]

And for small r (small u), the equation become,

s = - V / 4πTt

The measured values of mound dissipation (s) can be plotted against t^-1 to obtain the value of transmissibility T.

Initial head of aquifer

H(t)

Water level at start of test

Aquifer

Dipper Water Meter

The method can be illustrated by the following example (adapted from Ferris and Knowles 1963):

The following data is collected from a slug test. The volume of the slug is 0.148 m³.

-S (cm)	7.9	7.6	6.1	5.2	4.9	4.6	4.3	3.7	3.4	3.0	2.8	2.4	2.1	1.8	1.5	1.2	0.9
1/t min^-1	0.80	0.75	0.67	0.52	0.46	0.44	0.41	0.36	0.33	0.30	0.27	0.23	0.21	0.18	0.15	0.12	0.08

Plotting the values on coordinate paper and selecting arbitrarily a point on the line with coordinates –s = 6.3 cm and 1/t = 0.6 min^-1’ then,

T = V / 4 π t (s) = (0.148) (0.6) / [4 π (6.3 / 100)] = 0.11 m² / min.

0.0 0.2 0.4 0.6 0.8 1.0

1/t, min^-1

Location of Hydro-geologic Boundaries:

The existence of hydro-geologic formations such as recharge and impermeable boundaries limit the continuity of aquifers and limit the use of non-equilibrium formulas. Fortunately most surface recharge boundaries near pumping wells are visible and can be located before incorporating their effect on pumping. Unfortunately, impermeable subsurface boundaries such as faults, dykes and other formations are not apparent. The location and orientation of these barriers need to be established before any attempt is made to define an aquifer and its boundaries. Ferris (1962) employed pumping test data to define such boundaries.

As mentioned earlier, in the absence of hydro-geologic barriers, the Cooper-Jacob approximation would yield a straight line when plotted on semi logarithmic paper. However, when the propagating cone of depression encounters hydro-geologic barriers, the rate of drawdown changes and, depending on the nature of the barriers, the slope of the line will change. Since the radius of influence is defined as the radius, at which drawdown s is zero, then,

log (2.246Tt / r²S) = 0

and for constant T and S,

r_i² / t_i = r_r² / t_r

where r_i indicates the distance from the observation well to the image well, and r_r indicates the distance from the observation well to the pumping well.

When the propagating cone of depression reaches the observation well, drawdown begins and s vs. log t appears as a transitional curve. The curve becomes a straight line when u becomes small, and the properties T and S can be determined from the line using Cooper-Jacob method as described earlier. As the pumping continues, the effect of the boundary will eventually be felt at the observation well due to the reflection of the cone of depression. From this point on, the observed drawdown is the net drawdown due to both the pumping well and the image well causing the straight line to become steeper.

Semi log t

IW

Location of Impervious Boundary:

1) Draw s vs. t on semi log paper.

2) Choose arbitrary drawdown from the first segment of s vs. t line at the region where there is no boundary effect.

3) For the same value of ∆s₁, find the time (t_i) it takes for an image well to produce the same ∆s₁ at the region where there is boundary effect.

4) Find the distance, r_i to the image well from the equation r_i² / t_i = r_r² / t_r.

5) The distance r_i defines a radius of a circle at which image well lies. In order to define the location uniquely, two more observation wells are needed and the location of the image well will be the point of intersection of the three arches (see figure below)

Flow near Coastal Aquifers

Seawater encroachment into coastal aquifers pause great problems in populated urban areas that border the sea such as Dammam, Al-Khubar and Jeddah. Generally, there is a fresh groundwater gradient causing flow toward the sea. This gradient is due to excess fresh water in the coastal aquifer. The existence of this gradient is balanced by heavier, underlying seawater, and under natural, undisturbed conditions, a state of equilibrium exists between the two water bodies. Under ideal conditions, a zone of contact separates the freshwater body from the heavier salt water that lies underneath. This zone, known as interface, is actually a transition zone caused by dispersion of salt water. The density across this zone varies from the density of fresh water to the density of seawater. When conditions at the vicinity of the coast change, due to pumping or recharge, the water table (or piezometric surface) decreases or increases and the seawater interface advances inland or retreats until a new equilibrium state is reached. The advancement of seawater lens due to access pumping in the coastal areas is referred to in groundwater literature as seawater intrusion. When pumping becomes excessive, the intruding salt-water zone widens resulting in upcoming towards the pumping wells. If the rate of pumping is not carefully controlled, the pumping wells will become contaminated with salt water, resulting in eventual abandonment of well field.

h_s

h_f

Ghyben-Herzberg Approximation:

Assuming stationary seawater and hydrostatic pressure distribution in the fresh water region, the hydrostatic balance between the two regions is illustrated by the U tube shown above. In fact, the situation is actually a hydrodynamic equilibrium rather than hydrostatic one because fresh water flows toward the sea. It can be shown, however, that the Ghyben-Herzberge approximation that is based on hydrostatic consideration gives satisfactory results for locating salt-fresh water interface as long as the slope of the water table remain small. This approximation, however, is not valid at areas near the shoreline where Dupuit assumption is not valid.

Location of Salt-Fresh Water Interface:

z

Interface

η

Salt water ρ_s

Fresh water

ρ_f

If the elevation head is z, the head at fresh water and salt water

are defined as,

h_f = P_f / γ_f + z and,

h_s = P_s / γ_s + z

Since the pressure is continuous across the interface, the pressure

Is the same at the interface (at elevation η), then,

ρ_f (h_f – η) = ρ_s (h_s – η)

The equation become,

η = ρ_s / ∆ρ h_s – ρ_f / ∆ρ h_f

where, ∆ρ = ρ_s - ρ_f

For the special case of static salt water, ∆η can be written as,

∆η = – ρ_f / ∆ρ ∆h_f

The above equation relates the change of interface elevation with the change of fresh water head at the interface. If Dupuit assumption applies, ∆ h_f ≈ ∆h where ∆h is the change of elevation of water table. If z is the distance below sea level, the equation become,

z = (ρ_f / ∆ρ) h

And for γ_s = 1.025 and γ_f = 1.000 gr/cm³ then, z = 40 h_f

The above is known as the Ghyben-Herzberge approximation, which states that the depth of the interface at any point below sea level is approximately 40 times the height of the fresh water table above it.

40∆h

h_f

h_s = 40 h_f

Equipotential lines

Interface

Actual depth of interface

The discharge Q of fresh water is defined by Darcy’s law,

Q = K (z + h) dh/dx

Using the equation Ghyben-Herzberg Approximation for z and integrating, subject to h = 0 @ x = 0, the equation become,

h = {2QΔρ / [K (ρ + Δρ)] x }^1/2 and for small ∆ρ,

h ≈ {[2QΔρ / (K ρ) ]x }^1/2

This is the equation for water table elevation relative to sea level. Using Ghyben-Herzberg Approximation again,

h = (∆ρ/ρ) z = {[2QΔρ / (K ρ) ]x }^1/2

z = ρ / Δρ {2QΔρ / [K (ρ + Δρ)] x }^1/2, and for: (ρ + ∆ρ ≈ ρ)

z ≈ {2Qρ / (K∆ ρ) x }^1/2

Which is the equation defining the interface (McWhorter et al 77). In the above analysis, the consequence of using Dupuit assumption, allows no fresh water outlet to the sea, a physical impossibility. As an alternative, Glover (1959) provided better approximation for the interface elevation,

z = {[2Qρ / (K∆ ρ)] x + (ρQ/∆ρK) ²}^1/2

and the equation for the outlet (@z = 0) is,

x₀ = - Qρ / 2K∆ ρ

EXAMPLE:

Water is discharged to the sea at a rate of 10^-5 m³ / sec. The impervious stratum below sea level is 50 meters. If the hydraulic conductivity K is 5 x 10^-5 m/s, a) Calculate the toe position of the interface, b) calculate the amount of salt-water intrusion if the fresh water flow is reduced by 50%. Assume ρ_s = 1.05

	L

From Glover equation, @ z = 50 m

L = {K∆ ρ / 2Qρ [z² – (ρQ/∆ρK) ²}

= {(5x10^-5)(0.025) / [(2) (10^-5) ((1)] {(50)² – [(1) (10^-5)] / [(0.025) (5x10^-5)]²}

= [0.0625] [2500 – 64] = 152.25 m

The size of the fresh water outlet is,

X₀ = - Qρ / 2K∆ ρ

= - 10^-5 (1) / (2x10^-4x0.025) = - 2 m (negative because it is to the left x-axis)

Using the approximate equation’

L = (K∆ ρ) / 2Qρ z² = [(5x10^-5)(0.025)] / [(2x10^-5)(1)] (2500) = 156.25 m

Reducing Q by 50%,

L = {(5x10^-5)(0.025) / [(2) (0.5x10^-5) ((1)] {(50)² – [(1) (0.5x10^-5)] / [(0.025) (5x10^-5)]²}

= [0.125] [2500 -16] = 310.5 m

One Dimensional Solution of Flow Problems:

The differential equation describing flow in one dimension is,

∂s /∂t = C ∂²s / ∂x²

Subject to appropriate boundary and initial conditions. Examples of flow problems that is defined by above equation are,

1) Bank storage.

2) Line source.

3) Flow between two parallel drains.

4) Flow through earth dams.

a) Bank Storage: (Prescribed Drawdown)

Bank storage problems are encountered whenever an aquifer is hydraulically connected to streams and reservoirs. Depending on the elevation of water table relative to the stream elevation, the stream may recharge into the aquifer or feed from it. During flood period, the groundwater level is temporarily raised and then released to the adjacent streams following the storm. The volume of water released into the stream is referred to as bank storage.

The one-dimensional equation describing the idealized flow shown above is,

∂²s / ∂x² = (S/T) ∂s/∂t

The boundary and initial conditions, s (x,0) = 0, s (∞,t) = 0, s (0,t) = s₀

Introducing the variable,

u = x/ [4(T/S)t]^1/2 the equation becomes an ordinary differential equation,

d²s/du² + 2u ds/du = 0

The modified boundary conditions are s = 0 @ u = 0 and s = 0 @ u= ∞

s = s₀ [ 1 – erf (u) ] where, u = x/ [4(T/S)t]^1/2

Which is the response of the aquifer to a unit step drawdown of stream stage at t=0. (Hall and Moench, 1972).

The discharge to the channel per unit length (assuming abrupt change in water stage) is,

Q = -T [∂s / ∂x]_x=0

= -s₀ d/dU (erf u ) [∂u/∂x]_x=0

This yield

Q = s₀ T/ (π (T/S) t]^1/2

It is important to note that the above analysis is based on ground water response to discrete change of stream stage. Analysis based on continuous change of stage requires the integration of ds/dτ from t=0 to t=t where τ is time step function, the integration yields (Cooper and Rorabaugh 1963),

s = ʃ_0-t (ds₀/dτ) erfc [x/√4(T/S)(t-τ)] dτ, and,

Q ={T/[π(T/S)]^1/2} ʃ_0-t (ds₀/dτ) (t-τ) ^-1/2 dτ

EXAMPLE

A canal adjacent to an alluvial aquifer with water table elevation the same the bottom of the canal. Water is diverted to the canal at the following rates,

a) A sudden channel drawdown of 0.5 meters.

b) Continuous channel drawdown with the following schedule: -s₀ = 0.05 t meter for 0< t < 60 minutes

If T = 5 m² / min and S = 0.05 Find the equation for discharge of aquifer into the channel vs. time.

a) For sudden channel drawdown of 3 meters,

s = s₀ {1 – erf(x / [4(T/S)t ]^1/2}

Q = - T ∂s/x l_x=0

Q = s₀ T / [4(T/S)t ]^1/2 = (3) (5) / [4x100xt]^1/2

Q = 0.11 t^1/2

The equation for the discharge,

Q =( T/√π(T/S)) ʃ_0-t (ds₀/dτ) (t-τ) ^-1/2 dτ

For 0 < t < 60:

Q = -(2)[5/(π(5/0.05)]^1/2ʃ_{0→ 60} (0.05 (t-τ) ^-1/2 (Q is multiplied by 2 to account for both sides of the channel)

Q = -(2)[(5/(π(5/0.05)]^1/2 (0.05) (2)(t- τ)^1/2 l_{0 to t}

Q = - (2)[(5/(π(5/0.05)]^1/2 (0.05) (2)(t)^1/2

Q = - 0.0564 t^1/2

60 m t

0 30

b) Bank Storage: (Prescribed Discharge)

The one-dimensional equation for Darcy’s velocity is,

∂²q / ∂x² = (S/T) ∂q/∂t

The boundary and initial conditions for Darcy’s velocity is,

q(x,0) = 0, q(∞,t) = 0 and q(∞,t) = Q/2b where b is the saturated thickness. The solution is,

q = Q/2b erf {1- x / [4(T/S)]^1/2]

Since q = K ∂s/∂x, substituting and integrating and using the term u = x / √4(T/S)t yield,

s = [Q x / (2T√π)] {exp (-u²) / u + erf (u) - 1}

And, at x = 0, the drawdown is,

s₀ = Q [π(T/S)t)^1/2 / (πT)

The above equations are plane-source equations. They can be applied to channel seepage into underlying aquifers.

EXAMPLE:

A horizontal water table exists 10 meters below a canal. When water is diverted into the canal, leakage occurs at an average of 0.0017 m³/minute per meter. A) Estimate the time at which the water table will have risen to 2 meters above the original level and B) the height of the mound 20 meters from the canal. The properties of the aquifer are: T=0.2 m³/min and S_y = 0.05.

The time at which the water table has built up beneath the stream with s₀ = -2 m is,

t = (πTs₀/Q)² [1/π(T/S)]

t = [πx0.2x(-2)/(-0.0017)]² (1/4π)

= 43482.3 min = 30.2 days

2m

The height of the mound 50 meters from the canal:

At x = 50 , u = x/√4(T/S) t

u = (50) /√4(4)(43482.3) = 0.06

s = [Qx / [(2) (T)(√π)] {exp (-(0.06)²) / 0.06 + erf (0.06) – 1}

= [0.0017 / (2) (0.2) (√π) {0.996 / 0.06 + 0.156 - 1} = 0.4

= 0.0315 m