Hillel, pp. 190 - 192

12

Predictive Models for Saturated Hydraulic Conductivity

CE/ENVE 320 – Vadose Zone Hydrology/Soil Physics Spring 2004 Copyright © Markus Tuller and Dani Or 20022002-2004

The Kozeny-Carman Predictive Model for Ksat ● Efforts to predict hydraulic conductivity based on soil structural attributes such as porosity, particle size, etc. date back to early 1900. ● The idea is extremely attractive because of the obvious connection between water flow and the distribution of solids and pore spaces in the soil medium. ● In 1927 Kozeny introduced a model relating soil permeability to porosity. Permeability

ρ g η

Fluidity

⎛ρ g⎞ K s = k ⎜⎜ ⎟⎟ liquid density [M L ] ⎝ η ⎠ -3

acceleration of gravity [L t-2] viscosity coefficient of the liquid in [Pa s] [M L-1 t-1]

Fluidity at 20oC is: 9.81x106 [m s]-1. Note that permeability has dimensions of L2. Copyright© Markus Tuller and Dani Or2002-2004

The Kozeny-Carman Predictive Model for Ksat • Kozeny’s model was based on application of Poiseuille’s law to flow through granular porous media lacking structure or consolidation. • He idealizes the soil porous system as a bundle of cylindrical capillaries with a single representative radius. • Rearranging Poiseuille’s law for laminar flow in tubes, considering

the fraction of pores given by porosity n, yields Kozeny’s model in its simplest form:

n⎛R⎞ k= ⎜ ⎟ 2⎝ 2 ⎠

2

2

[m ]

A R2 π h RH = = P 2 Rπ h RH= Hydraulic Radius (ratio pore volume to wetted surface area) n Porosity Copyright© Markus Tuller and Dani Or2002-2004

The Kozeny-Carman Predictive Model for Ksat Kozeny’s approach was extended by Carman to take the form: 3

n k= 2 2 cτ Av (1 − n) 2

2

[m ]

c Empirical coefficient representing particle size and shape (2-3 for non-cubic particles) τ Tortuosity of the flow path (21/2 for main flow direction 45o from the vertical axis) Av Specific surface per unit volume. Above equation can be simplified by introducing a hydraulic radius RH:

RH = Av (1 − n )

n3 k= 2 2 cτ RH Copyright© Markus Tuller and Dani Or2002-2004

The Kozeny-Carman Predictive Model for Ksat Tortuosity: Ratio of the average roundabout path of a fluid molecule to the apparent or straight flow path. Tortuosity is a dimensionless geometric parameter of porous materials, which, though difficult to measure precisely, is always greater than 1 and may exceed 2.

Copyright© Markus Tuller and Dani Or2002-2004

The Kozeny-Carman Predictive Model for Ksat k=

n3 cτ 2 R h 2

The product cτ2 is commonly about 5. The pore shapes, surfaces, and tortuosity are interdependent, thus the product cτ2 Rh2 is often lumped into a single “formation factor”, F2.

n3 k= 2 F Ks is estimated as:

⎛ρg⎞ n ⎛ρg⎞ ⎟⎟ = 2 ⎜⎜ ⎟⎟ K s = k⎜⎜ ⎝ η ⎠ F ⎝ η ⎠ 3

Copyright© Markus Tuller and Dani Or2002-2004

The Kozeny-Carman Predictive Model for Ksat The Kozeny-Carman approach provides unsatisfactory estimates of Ks in many soils, due to the assumption of uniform pore radii. Works well for sands and other materials with uniform pore size distribution.

Copyright© Markus Tuller and Dani Or2002-2004

Revil & Cathles Predictive Model for Ksat Revil and Cathles [1999] developed a predictive model for permeability using concepts from electrical conductance in porous media. They use the concept of a formation factor F that is a dimensionless, scale invariant parameter that characterizes the pore space topology. F is related to the porosity n and pore geometry by the empirical Archie relationship:

F=n

−m

m is the so-called “cementation exponent” and varies with pore geometry in a range of 1 to 4. (ratio of pure solution electrical conductivity to the bulk electrical conductivity of the saturated porous medium.

Copyright© Markus Tuller and Dani Or2002-2004

Revil & Cathles Predictive Model for Ksat Cementation Exponent m For media with perfectly spherical grains m=1.5 For sands with porosities between 0.03 and 0.3 m=1.5 to 2.0 For media where the pore space consists entirely of open interconnected fractures and cracks m=1.1 to 1.3 (most of the porosity conducts electrical current) If large pores are connected to narrow throats m>2.5

The derivations of Revil and Cathles [1999] yielded a simple expression:

d 2 n 3m k= 24 d is the mean grain diameter (from particle size analysis) of granular materials. Copyright© Markus Tuller and Dani Or2002-2004

Revil & Cathles Predictive Model for Ksat The approach yields remarkable match for granular materials, and was later extended to predict hydraulic conductivities of sand-clay mixtures.

1 mD (milidarcy) = 1x10-15 m2 ~10-8 m/s (water 200 C) Copyright© Markus Tuller and Dani Or2002-2004

Characterization of Flow Conditions Flow regimes in a porous medium are characterized according to: 1) Changes in flow attributes in time and space (steady vs. non steady flows) 2) Saturation state of the system (saturated vs. unsaturated) 3) Geometry of flow (one-dimensional or multi-dimensional). Constant Head Method for Ksat: Steady-state, saturated, and one-dimensional flow conditions. (water flux and other flow attributes such as hydraulic potentials and water content do not change over time at any given location).

Falling Head Method for Ksat: Non steady, saturated, and one-dimensional flow conditions. (Flux density and hydraulic potentials vary with time.) Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Flow

CE/ENVE 320 – Vadose Zone Hydrology/Soil Physics Spring 2004 Copyright © Markus Tuller and Dani Or 20022002-2004

Flow of Water in Unsaturated Soils • Most flow processes in the field occur under unsaturated conditions, or in the presence of an air-phase. • This drastically modifies water flow pathways, with decreasing water content flow through water filled pores become constrained to narrower and tortuous channels. • Darcy's law, originally conceived for saturated flow only, was extended by Edgar Buckingham (1907) to unsaturated water flow. Buckingham's main assumption was that unsaturated hydraulic conductivity is a function of water content or matric potential, i.e. the degree of wetness.

⎡ ∆H ⎤ ⎡ ∆(h + z ) ⎤ ⎡ ∆h ⎤ Jw = −K(h)⎢ = −K(h)⎢ = −K(h)⎢ + 1⎥ ⎥ ⎥ ⎣ ∆z ⎦ ⎣ ∆z ⎦ ⎣ ∆z ⎦

⎡ dh ⎤ Jw = −K(h) ⎢ + 1⎥ ⎣ dz ⎦ Unsaturated Hydraulic Conductivity expressed as a function of Matric Potential h

Copyright© Markus Tuller and Dani Or2002-2004

Saturated versus Unsaturated Flow The driving force for flow in a saturated soil under constant gravitational potential is the gradient of positive pressures, whereas under unsaturated conditions the gradient is of negative matric potentials. When the soil is completely water-saturated, the pores conduct at a maximal rate Ks. When the soil desaturates, some of the water-filled pores become partially or completely air-filled, and the cross-sectional area available for mass flow decreases K(h).

Equations for saturated flow have a constant saturated hydraulic conductivity (Ks) Unsaturated flow equations have an unsaturated hydraulic conductivity function, K(h).

Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Hydraulic Conductivity Function crossover

● The unsaturated hydraulic conductivity is a highly nonlinear function of water content, K(θ), or matric potential, K(h)often spanning more than 8 orders of magnitude! ● At saturation a sandy soil has higher conductivity than clay due to larger water-filled pore spaces. ● At lower suctions large pores in sandy soil drain resulting in a steep decrease in conductivity. ● Loamy and clayey soils, having smaller (on average) pores, retain greater number of water-filled pores resulting in higher hydraulic conductivity than sandy soils at low potentials.

Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Hydraulic Conductivity Function Estimates for unsaturated hydraulic conductivity for various soil textural classes can be obtained from the UNSODA and NRCS Soil Survey databases. 4

4 2

1

a

2

2

3

2

12

11 8

0

1

b 7

9

0

10

8

6

5

9

12

-2

-2

3 4

-4

lo g K (cm /d )

log K (cm /d)

11

-6 -8

-4 -6 -8

UNSODA -10 -12

6

-14 7

-16 0.0

1 Sand 7 Sandy clay loam 2 Loamy sand 8 Clay loam 3 Sandy loam 9 Silty clay loam 4 Loam 11 Silty clay 5 Silt 12 Clay 6 Silt loam

4

0.1

5

0.2

0.3 θ

0.4

0.5

Soil Survey

-10

1 Sand 7 Sandy clay loam 2 Loamy sand 8 Clay loam 3 Sandy loam 9 Silty clay loam 4 Loam 10 Sandy clay 5 Silt 11 Silty clay 6 Silt loam 12 Clay

-12 -14 -16 0.0

0.1

0.2

0.3

0.4

0.5

θ Copyright© Markus Tuller and Dani Or2002-2004

Parametric Models for K(θ) Some researchers have extended the concepts discussed in the previous section and derived closed-form analytical predictions of the unsaturated hydraulic conductivity. The basic step is replacing the summations with integrals of SWC analytical (parametric) expressions The van Genuchten and Brooks and Corey SWC Models are amenable to such integration. Van Genuchten (1980) derived an expression for the unsaturated hydraulic conductivity based on the same parameters he used for the SWC curve and on a modification of Jackson’s model proposed by Mualem (1976). The resulting water content-dependent unsaturated hydraulic conductivity K(θ) is given by:

⎡ ⎛ 1 ⎞ K (Θ ) = Θ ⎢1 − ⎜1 − Θ m ⎟ ⎢ ⎜ ⎟ Ks ⎠ ⎢⎣ ⎝ 1 2

m ⎤2

⎥ ⎥ ⎥⎦

θ − θr Θ= θs − θ r Copyright© Markus Tuller and Dani Or2002-2004

Parametric Models for K(h) van Genuchten-Mualem (VGM) model: The unsaturated hydraulic conductivity can also be expressed in terms of matric head (h):

K (h ) = Ks

[

⎡⎛ 1 − (αh ) n −1 1 + (αh ) ⎢⎣⎜⎝

[1 +

]

n −m

⎞⎟⎤ ⎠⎥⎦

2

]

m (αh ) n 2

Copyright© Markus Tuller and Dani Or2002-2004

Parametric Models for K(h) Brooks and Corey (BC) model: Similarly expressions can be derived for the Brook's and Corey model.

K (h ) ⎛ h = ⎜⎜ Ks ⎝ hb

⎞ ⎟⎟ ⎠

−2−3λ

where hb is the matric head at air-entry value or "bubbling pressure", and λ is the pore distribution index parameter. Expressed as K(θ):

K (Θ ) =Θ Ks

3+

2

λ

θ − θr Θ= θs − θ r

Copyright© Markus Tuller and Dani Or2002-2004

Parametric Models for K(θ) & K(h) As input information for the VG and BC models for unsaturated hydraulic conductivity we need to measure a few h-θ pairs (SWC) and saturated hydraulic conductivity Ks: (1) We fit the parametric expressions for the SWC to measured data to obtain free model parameters α, n, θr (VG) or hb, λ (BC) (2) We use measured Ks and the free model parameters from the SWC to predict unsaturated hydraulic conductivity

Copyright© Markus Tuller and Dani Or2002-2004

Parametric Models for K(h) Another widely used parametric model for unsaturated hydraulic conductivity was proposed by Gardner (1958):

K (h ) = K s e

βh

Where β (dimensions of L-1) is a parameter related to typical pore size of the medium. The primary advantage of Gardner’s model is that it allows analytical solutions to certain flow problems.

Why do we need parametric models? •

For obtaining analytical solutions to flow and transport problems

•

For comparison between different types of soils or other porous media. Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow in Soils

CE/ENVE 320 – Vadose Zone Hydrology/Soil Physics Spring 2004 Copyright © Markus Tuller and Dani Or 20022002-2004

Steady State Unsaturated Flow Steady-state unsaturated flow may develop under a variety of conditions. (1) When a constant flux at the soil surface, less than the soil's Ks is maintained for extended periods of time by use of a sprinkler or other means. (2) The evaporation rate from a shallow water table may also attain (approximately) a steady state. For vertical steady-state unsaturated flow we may write the Buckingham-Darcy equation:

⎡ dh ⎤ J w = − K (h ) ⎢ + 1⎥ ⎣ dz ⎦

dh = −dz Jw 1+ K (h )

This ordinary differential equation may be integrated to yield the integral form of the Buckingham-Darcy equation (next slide). Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow Integral form of the Buckingham-Darcy equation:

dh = −dz J 1+ w K (h ) h2

z2

dh ∫ J w = − ∫ dz = z1 − z 2 h1 1 + z1 K (h ) E.g., Gardner’s Model K (h ) = K s e β h

When K(h) has a known functional form, the integral on the left side may be evaluated exactly to produce an analytical solution for the water flux or for the distribution of the matric potential in the soil profile, h(z). Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow Analytical Solution for Unsaturated Steady State Flow to a Water Table: Given a constant flux (Jw < Ks) from the soil surface to a water table at a depth L. We assume that the unsaturated hydraulic conductivity, K(h), is given by Gardner's model K(h)=Ks eβh We want to find the distribution of soil matric potential in the soil profile from the soil surface to the water table. (1) The first step is to introduce the parametric expression of K(h) into the integral form of the Buckingham-Darcy law

h2

z

2 dh ∫ J w = − ∫ dz = z1 − z 2 h1 1 + z1 K (h )

h (−L)

∫

h (z)

dh J w −β h 1+ e Ks

−L

= − ∫ dz = L + z z

Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow (2) We set z=0 at the soil surface, and z=-L at the water table (where h=0). From a table of integrals we find the following solution:

[

]

dx x 1 ax = + ln p + qe +C ∫ p + qeax p ap

(3) Applying this solution to our problem with the appropriate boundary conditions yields:

⎡ J w −β h (z) ⎤ 1 ⎡ J w ⎤ z + L = h (z) − ln ⎢1 + e ⎥ + ln ⎢1 + ⎥ β ⎣ Ks β K s⎦ ⎦ ⎣ 1

(4) Rearranging and log-transforming yields the desired form in terms of h(z):

h ( z) =

⎡ J ⎤ J ln ⎢− w + e − β ( L+ z ) + w e − β ( L + z ) ⎥ Ks β ⎣ Ks ⎦ 1

One can verify that at z=-L, h(z=-L)=0. The term in the square brackets is limited to values between 0 and 1. This is important for analyses of maximal upward flow from a water table. Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow Graph showing matric potential distributions within the soil profile for various potentials (previous equation).

Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow - Example Problem Statement: Compute and plot the vertical distributions of the soil matric potential under steady flow into a water table at a depth of 3 m from the surface for the flow rates of Jw/Ks= -0.001, -0.01, -0.1, and -0.5. Note that the sign of the downward flux is negative. The unsaturated hydraulic conductivity of the soil is based on Gardner's (1958) model, and is characterized by β=5 m-1. Estimate the hydraulic gradient (dH/dz) from 0 to 1m below the soil surface. Solution: 1. The depth dependent matric potential may be calculated according to the following relationship:

⎡ Jw J w −β (L+z) ⎤ − β ( L+ z) +e + h (z) = ln ⎢− e ⎥ Ks β ⎣ Ks ⎦ 1

Results are listed on the next slide. Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow - Example Tabulated Results:

⎡ Jw J w −β (L+z) ⎤ − β ( L+ z) +e + e h (z) = ln ⎢− ⎥ Ks β ⎣ Ks ⎦ 1

Jw/Ks Z [m] 0.00 -0.25 -0.50 -0.75 -1.00 -1.25 -1.50 -1.75 -2.00 -2.25 -2.50 -2.75 -3.00

-0.001 h(z) [m] -1.3815 -1.3813 -1.3808 -1.3790 -1.3727 -1.3522 -1.2936 -1.1667 -0.9725 -0.7419 -0.4978 -0.2495 0.0000

-0.01 h(z) [m] -0.9210 -0.9210 -0.9210 -0.9208 -0.9201 -0.9179 -0.9104 -0.8861 -0.8188 -0.6805 -0.4788 -0.2451 0.0000

-0.1 h(z) [m] -0.4605 -0.4605 -0.4605 -0.4605 -0.4604 -0.4602 -0.4595 -0.4571 -0.4487 -0.4221 -0.3499 -0.2055 0.0000

-0.5 h(z) [m] -0.1386 -0.1386 -0.1386 -0.1386 -0.1386 -0.1386 -0.1385 -0.1382 -0.1373 -0.1340 -0.1229 -0.0882 0.0000 Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow - Example Solution – Graph:

- 1 .6

- 1 .4

- 1 .2

S o il M a t r ic P o t e n t ia l ( m ) -1 - 0 .8 - 0 .6

- 0 .4

- 0 .2

0

0 .0 0

0 .5 0

Soil Depth (m )

1 .0 0

1 .5 0

2 .0 0

2 .5 0

3 .0 0

3 .5 0

Copyright© Markus Tuller and Dani Or2002-2004

Steady State Unsaturated Flow - Example Solution - Continued: 2. To estimate the hydraulic gradient we have to solve:

dH d (h + z ) dh = = +1 dz dz dz J w K s = −0.001

dH − 1.38149 + 1.37268 + 1 = 0.99119 = dz 1

J w K s = −0.01 ⇒ dH dz = 0.99911 J w K s = −0.1 ⇒ dH dz = 0.99992 J w K s = −0.5 ⇒ dH dz = 0.99999 Copyright© Markus Tuller and Dani Or2002-2004

Prediction of Unsaturated Hydraulic Conductivity The relationships between the hydraulic conductivity and the degree of saturation or matric potential arouse interest in methods for predicting the soil hydraulic conductivity based on SWC. Most models are based on simplifying assumptions: 1) The complicated soil pore geometry is represented by a bundle of continuous cylindrical capillaries. 2) The flow is assumed to occur only in these capillaries (no bypass flow) 3) The length of the bundle of capillaries (Lc), is assumed to be larger than the length (L) of the soil column it represents; the ratio Lc/L is the tortuosity factor (τ). Soil Pores

Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity Poiseuille's law for laminar flow in cylindrical tubes is applied to the bundle of capillaries with a gradient ∆H across its length. This enables computation of the total volume of flow through the column by: Soil Pores

Poiseuille's Law

π R 4 ∆P Q= 8η L

π ρ w g ∆H M 4 = ∑ N jQ j = N R ∑ j j η 8 L c j=1 j=1 M

Q tot Nj

is the number of capillaries of radius Rj in the bundle

M

is the number of different (arbitrarily assigned) capillary size classes Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity The flux (Jw) through the column is:

Jw nj

Q tot π ρ w g ∆H M 4 = = n R ∑ j j A 8η L c j=1

is the number of tubes with radius Rj per unit area (nj =Nj/A)

The next step is to convert the soil water characteristic curve into a soil water content vs. equivalent pore radii curve. This is done by using the capillary-rise equation to convert a matric potential value to an equivalent capillary radius (see next slide).

Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity 1) The total saturation water contend θs is divided into a certain number M of equal increments ∆θv where the matric potential hj corresponds to a water content of θs- jx∆θv 2) We assume that all tubes with radii greater than Rj have drained at matric potential hj

− 2σ r>Rj = ρw g h j

Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity SWC-Curve

1.E+04

1.E+02 1.E+01 1.E+00 1.E-01 1.E-02 0

0.1

0.2

0.3

0.4 3

Water Content- Equivalent Capillary Radius Curve

0.5

3

Volumetric Water Content [m /m ]

1.E-02 1.E-03

Capillary Rise Equation

− 2σ Rj = ρw g h j

Equivalent Radius [m]

Matric Potential [-m]

1.E+03

Empty

1.E-04 1.E-05 1.E-06 1.E-07

Full

1.E-08 1.E-09 0

0.1

0.2

0.3

0.4 3

0.5

3

Volumetric Water Content [m /m ]

Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity

3) The number of capillaries having a radius Rj per unit area in each water content interval is: nj=∆θv/πRj2 (where ∆θv is interpreted as the fraction of the water-filled cross sectional area that is reduced when all capillaries having a radius of Rj drain).

Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity 4) Now we can plug in the number of capillaries nj into Poiseuille’s law and introduced a multiplication by -L/∆z=1 (where ∆z=0-L) to obtain a form similar to Darcy's law : 2 ⎡ π ρ w g ∆H M ∆θ v σ ∆θ v M 1 ⎤ L ∆H 4 ⎥ R j = −⎢ Jw = ∑ ∑ 2 2 8η L c j=1 π R j ⎢⎣ 2η ρ w g j=1 h j ⎥⎦ L c ∆z

Saturated Hydraulic Conductivity

τ Tortuosity Factor

τ σ 2 ∆θ v M 1 Ks = ∑ 2ηρ w g j=1 h j 2 Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity As water content decreases from θs to θs- i x ∆θv the tubes larger than Ri drain and do not contribute to water flow. The effect on the unsaturated hydraulic conductivity, which is a function of the water content, is given by

τ σ 2 ∆θ v M 1 Ks = ∑ 2ηρ w g j=1 h j 2 τ σ 2 ∆θ v K (θ s − i∆θ v ) = 2ηρ w g

M

1

j=i +1

j

∑h2

This equation, originally Developed by Childs and Collis-George (1950), serves as the basis for many predictive models of unsaturated hydraulic conductivity, most of which are based on SWC information and assume a capillary bundle geometry. Copyright© Markus Tuller and Dani Or2002-2004

Predictive Models for Soil Hydraulic Conductivity Even with a detailed SWC, there is still an unknown, i.e. the tortuosity (τ). One way to circumvent this problem is to measure the saturated hydraulic conductivity (which is relatively easy to measure) and use it to solve for τ, assuming it is constant. One such formulation was proposed by Jackson (1972):

∆θ ⎛ ⎜θj + 2 K(θi ) = K s ⎜ ⎜ θs ⎜ ⎝

⎞ ⎟ ⎟ ⎟ ⎟ ⎠

c

M

∑

2 j + 1 − 2i

j =i M

∑

j =1

h j2 2j − 1 h j2

where θj=(θs-j x ∆θv), c is a constant usually taken as c=1, and the water content is divided to M equal increments (of θv)

Copyright© Markus Tuller and Dani Or2002-2004

Jackson Model - Example The matric potential-soil water content function (SWC) for Sarpy loam soil is given by the VG model expressed for the matric potential (h) as:

Θ=

1 ⎤n

1 ⎡ − 1 h = ⎢Θ m − 1⎥ α⎢ ⎥⎦ ⎣

θ − θr θs − θ r

θs = 0.4 m3 m3 θ r = 0.03 m3 m3 α = 2.79 m −1 n = 1.6 m = 1 − 1 n = 0.375 K s = 50 mm hr c =1 M = 10

Copyright© Markus Tuller and Dani Or2002-2004

Jackson Model - Example First we split the SWC into 10 equal increments and calculate water content θ and corresponding matric potential h for the midpoint of each increment. ⎡ 1 ⎢⎛ hj = ⎜1 − ⎟ α ⎢⎝ 2M ⎠ ⎢⎣

1 − 2 j −1⎞ m

1.E+04

Matric Potential [-m]

1.E+03

1 ⎤n

− 1⎥ ⎥ ⎥⎦

1.E+02 1.E+01 1.E+00 1.E-01 1.E-02 0

0.1

0.2

0.3

0.4 3

0.5

3

Volumetric Water Content [m /m ]

θ j = θs − (θs − θ r )

(2 j − 1) 2M

Copyright© Markus Tuller and Dani Or2002-2004

Jackson Model - Example Now we can apply Jackson’s equation and calculate K(θ) for the midpoints of each increment. If we consider for example the calculation of K(0.289):

K(0.289 ) = 50

*

0.289 13.0457 = 2.2038 0.4 213.844

[mm hr ]

j, i

hj (m)

1/hj2

θj+∆θ/2

Σ(2j+1-2i)hj2

K(θj)

1

0.1079

85.83

0.400

213.844

50.000

2

0.2446

16.72

0.363

70.8035

15.024

3

0.3919

6.51

0.326

30.3092

5.7757

4

0.5790

2.98

0.289

13.0457

2.2038

5

0.8424

1.41

0.252

5.27549

0.7771

6

1.253

0.637

0.215

1.89699

0.2384

7

1.982

0.254

0.178

0.56445

0.0587

8

3.556

0.079

0.141

0.12307

0.0101

9

8.430

0.0141

0.104

0.01515

9.21E-4

10

52.806

3.59E-4

0.067

0.00036

1.40E-5

j is a pointer whereas i is a counter

** ∆θ/2=0.37/(10x2)=0.0185; ∆θ/2=0.37/(10x2)=0.0185; note that hj corresponds to θj, which is less than the water content indicated in the table by the amount of ∆θ/2. ∆θ/2. Copyright© Markus Tuller and Dani Or2002-2004

Jackson Model - Example When we compare our results with the analytically derived van Genuchten-Mualem (VGM) K(θ) model (see next slides) we see good agreement.

Copyright© Markus Tuller and Dani Or2002-2004

Hillel, pp. 210 - 219

15

Unsteady and Unsaturated Flow – The Richards Equation

CE/ENVE 320 – Vadose Zone Hydrology/Soil Physics Spring 2004 Copyright © Markus Tuller and Dani Or 20022002-2004

Unsaturated Nonsteady Flow Unsaturated nonsteady flow is the most common flow regime under field situations, where water content θ, the matric potential h, and the flux Jw vary in time and space. The application of Buckingham-Darcy's law to describe the flow is no longer sufficient as the flux varies from one location to another, a condition called flux divergence. The difference between the flux entering a volume element of soil and that exiting the same volume element is the change in water content, or storage (CONSERVATION OF MASS = CONTINUITY PRINCIPLE).

Copyright© Markus Tuller and Dani Or2002-2004

The Richards Equation The principal of continuity may be formalized mathematically for one-dimensional vertical flow as:

∂θ v ∂J w = ∂t ∂z Note that the partial derivatives indicate possible variations of flux and water content in space and in time. The three-dimensional form of the continuity equation is given as:

⎛ ∂J wz ∂J wy ∂J wx ∂θ + = ∇ ⋅ J w = ⎜⎜ + ∂y ∂x ∂t ⎝ ∂z

⎞ ⎟⎟ ⎠

del operator (gradient of J) Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Nonsteady Flow – Richards Equation To obtain a flow equation for unsaturated nonsteady flow we have to combine Buckingham-Darcy's law for the flux (Jw) with the continuity equation. For one-dimensional vertical flow we obtain:

∂θ ∂ ⎡ ⎛ ∂h ⎞⎤ = ⎢K (h )⎜ + 1⎟⎥ ∂t ∂z ⎣ ⎝ ∂z ⎠⎦ For three-dimensional flow we obtain:

∂h ⎞ ∂K (h ) ∂h ⎞ ∂ ⎛ ∂θ ∂ ⎛ ∂h ⎞ ∂ ⎛ = ∇ ⋅ [K (h ) ∇H ] = ⎜ K (h ) ⎟ + ⎜⎜ K (h ) ⎟⎟ + ⎜ K (h ) ⎟ + ∂z ∂z ⎠ ∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ This equation is known as the RICHARDS EQUATION, and is the most general unsaturated flow equation used in practice. Note that this equation is based on the assumption of isothermal conditions (uniform temperature) and flow through a rigid matrix (no swelling or shrinking of the soil). Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated & Unsteady Horizontal Flow Unsaturated and nonsteady horizontal flow in one dimension may be represented by the simplest form of the Richards equation, with no gravitational component, as:

∂θ ∂ ⎡ ∂h ⎤ = K (h ) ⎥ ⎢ ∂x ⎦ ∂t ∂ x ⎣ In this case we have two variables, θ and h, which vary in time and space, but only one equation. In some cases it is possible to reduce the equation to only one variable by using the interrelationship θ−h (the SWC curve). For water infiltration into a homogeneous soil having uniform initial water content it is advantageous to eliminate the matric potential (h).

Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Nonsteady Horizontal Flow We first apply the chain rule to the matric potential gradient to obtain: ∂θ ∂ ⎡ ∂h ⎤ ∂h dh ∂θ K ( h ) = = ∂t ∂ x ⎢⎣ ∂x ⎥⎦

∂x

dθ ∂x

Reciprocal of the SPECIFIC WATER CAPACITY Cw(θ)

dθ C w (θ ) = dh Cw(θ) is the slope of the SWC curve at any particular water content.

Next we express K(h[θ]) = K(θ), and define a new hydraulic function, the Soil Water Diffusivity D(θ) as:

dh K (θ ) D(θ ) = K (θ ) = dθ C w (θ )

The Soil Water Diffusivity is the ratio of unsaturated hydraulic conductivity to the specific water capacity. Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Nonsteady Horizontal Flow Soil Water Diffusivity D(θ):

dh K (θ ) D(θ ) = K (θ ) = dθ C w (θ )

Since the K(θ) and Cw(θ) functions are highly nonlinear, the diffusivity is also nonlinear. However, the range of variation in D(θ) is much smaller than the range of K(θ) variation because the nonlinearity in numerator and denominator cancel out to some degree. The resulting water content or diffusivity form of the Richards equation for one dimensional horizontal flow is then:

∂ ⎡ ∂θ ⎤ ∂θ = D(θ ) ⎥ ⎢ ∂t ∂ x ⎣ ∂x ⎦

The driving force is the gradient in water content “water concentration”

This is not in accordance with the basic physical principle of flow occurring due to potential gradients !! Special care is required when applying the diffusivity form of the Richards equation!! Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Nonsteady Horizontal Flow

∂ ⎡ ∂h ⎤ ∂θ = K (h ) ⎥ ⎢ ∂t ∂ x ⎣ ∂x ⎦

∂ ⎡ ∂θ ⎤ ∂θ = D(θ ) ⎥ ⎢ ∂t ∂ x ⎣ ∂x ⎦ Copyright© Markus Tuller and Dani Or2002-2004

Unsaturated Nonsteady Horizontal Flow Graphical presentation of the K(θ), Cw(θ), and D(θ) functions CW(θ)

1.E+00

0.30 -1

]

1.E-01 1.E-02

Water Capacity [m

Hydraulic Conductivity [m /hr]

K(θ) 1.E-03 1.E-04 1.E-05 1.E-06 1.E-07

0.25 0.20 0.15 0.10 0.05 0.00

1.E-08 0

0.1

0.2

0.3

0.4 3

0.5

0

0.1

0.2

0.3

0.4 3

3

0.5

3

Volumetric Water Content [m /m ]

Volumetric Water Content [m /m ]

D(θ) 1.E+00 Diffusivity [m2/hr]

1.E-01

D(θ ) = K (θ )

1.E-02 1.E-03

dh K (θ ) = dθ C w (θ )

1.E-04 1.E-05 1.E-06 0

0.1

0.2

0.3

0.4 3

0.5

3

Volumetric Water Content [m /m ] Copyright© Markus Tuller and Dani Or2002-2004

Boltzmann transformation for horizontal flow

λ = xt

−1

2

∂θ dθ ∂λ λt −1 dθ = = ∂t dλ ∂t 2 dλ dθ ∂λ ⎤ d ⎡ dθ ⎤ ∂ ⎡ ∂θ ⎤ ∂ ⎡ D ( θ ) D ( θ ) D ( θ ) = = dλ ∂x ⎥⎦ d λ ⎢⎣ dλ ⎥⎦ ∂ x ⎢⎣ ∂x ⎥⎦ ∂ x ⎢⎣ Copyright© Markus Tuller and Dani Or2002-2004

ODE for Unsaturated Nonsteady Horizontal Flow

λ dθ d ⎡ dθ ⎤ D(θ ) ⎥ − = ⎢ dλ ⎦ 2 dλ d λ ⎣

Analytical solution with constant diffusivity D(θ)=D:

⎡

x ⎤ ⎥ ⎣ 2 Dt ⎦

θ ( x, t ) = θ i + (θ 0 − θ i ) erfc ⎢

Copyright© Markus Tuller and Dani Or2002-2004

Diffusivity Form of the Richards Equation Special care is required when applying the diffusivity form of the Richards equation!! It should be applied only to homogeneous soils, where there is a unique relation between h and θ, neglecting hysteresis, and to describe horizontal flow where there is no gravitational component. What happens in layered soils? We cannot apply diffusivity form of Richards equation!! (the derivative ∂θ/∂x is not defined across the discontinuity) Matric potential is continuous across layers

Water Content is discontinuous across layers Copyright© Markus Tuller and Dani Or2002-2004

Diffusivity Form of the Richards Equation The diffusivity form of Richards equation is applicable: In homogenous soils for conditions where the vertical gravitational component is very small in comparison with matric influences (translated into a water content "gradient"). Examples: 1) Early phases of infiltration into relatively dry soil, 2) Evaporation of water from the soil surface. For special cases the diffusivity may be assumed as constant over a small range of water content. This facilitates analytical solutions of the Richards equation.

Copyright© Markus Tuller and Dani Or2002-2004

Parametric Forms of the Diffusivity Function Diffusivity is defined as:

D(θ ) = K (θ )

dh dθ

We can express K(θ) with the van Genuchten – Mualem relationship and substitute with the VG model for Θ: 1 ⎡ 1 ⎞ ⎛ K (Θ ) ⎢ ⎜ 2 m = Θ 1− 1−Θ ⎟ ⎢ ⎜ ⎟ Ks ⎠ ⎢⎣ ⎝

m ⎤2

⎥ ⎥ ⎥⎦

(1 − m) K s D(Θ ) = Θ α m (θ s − θ r )

⎤ θ − θr ⎡ 1 Θ= =⎢ ⎥ θs − θ r ⎢⎣1 + (α h ) n ⎥⎦

1 1 − 2 m

m

⎡ ⎢(Θ ) h=⎣

1 − m

α

1 ⎤n

− 1⎥ ⎦

m ⎡⎛ ⎤ 1 ⎞−m 1 ⎛ ⎞ ⎢⎜1 − Θ m ⎟ + ⎜1 − Θ m ⎟ − 2⎥ ⎢⎜ ⎥ ⎟ ⎜ ⎟ ⎠ ⎝ ⎠ ⎢⎣⎝ ⎥⎦

Copyright© Markus Tuller and Dani Or2002-2004

Parametric Forms of the Diffusivity Function Applying the same approach for the Brooks and Corey model yields:

Ks h b D(Θ ) = Θ λ (θ s − θ r )

2+

1

λ

Copyright© Markus Tuller and Dani Or2002-2004

Numerical Solutions of Richards Equation ● The high degree of nonlinearity in the coefficients of the Richards equation for unsaturated flow makes it impractical to obtain analytical solutions. ● It might be considered as one of the most challenging equations in environmental physics (considerably different than flow of heat or electricity). ● Except for a few special cases, obtaining solutions to practical flow problems represented by the Richards equation is feasible by means of successive numerical approximations, generally using computers. Copyright© Markus Tuller and Dani Or2002-2004

Numerical Solutions of Richards Equation Numerical Approximations (Finite Difference Method) are based on: 1) Discretization of the flow domain into a grid: "layers" for one dimensional flow or a mesh for two-dimensional flow. 2) Discretization of the continuous partial differential equation into a series of approximated difference equations, one for each of the K spatial nodes.

i −1 / 2

θ i j +1 − θ i j ( hi j − hi −j 1 + ∆z ) K i +1/ 2 ( hi +j 1 − hi j + ∆z ) − ∆z ∆z ∆z

∆t

3) Application of initial and boundary conditions. 4) Formulation of the problem into a computer code which is able to solve the equations for the entire flow domain at various time steps. Copyright© Markus Tuller and Dani Or2002-2004

Numerical Solutions of Richards Equation Example: A Numerical Approximation to the Richards Equation Find a numerical approximation to changes in water content with time under vertical unsaturated and nonsteady flow into a soil with a known initial water content, and known hydraulic functions θ(ψ) and K(ψ). Solution: We first discretize the soil profile (the spatial domain) into units of ∆z, and the time into units of ∆t. We denote the water content at the i-th depth increment and at the j-th time step as θij; similarly we express the corresponding matric head as hij (see graph previous slide) The hydraulic conductivity, which is a function of the matric head (or of the water content), relates transfer between two soil layers and thus must be averaged to represent both layers as: K([hij + hji+1]/2) = Kji+1/2. The resulting numerical approximation for BuckinghamDarcy’s law for vertical flow is then given by: j j j j j j j+1 j ∂θ θ i − θ i ∂ ⎡ ⎛ ∂h ⎞⎤ (h i −1 − h i + ∆z) K i−1 2 (h i − h i+1 + ∆z) K i+1 2 − = = ⎢K ( h ) ⎜ + 1⎟⎥ = 2 ∆t ∂t ∂z ⎣ ∆z ∆z 2 ⎝ ∂z ⎠⎦ Copyright© Markus Tuller and Dani Or2002-2004

Numerical Solutions of Richards Equation Example – continued: Combining terms and rearranging allows us to solve for the only unknown in this equation, which is θij+1, the water content for the i-th depth increment at the next (future) time step:

θ ij+1 = θ ij +

[

∆t j j j j j j ( h − h + ∆ z ) K − ( h − h + ∆ z ) K i −1 i i i +1 i −1 2 i +1 2 2 ∆z

]

Because this is a discrete approximation of a continuous process, increments ∆z and ∆t should be kept small. In addition, initial soil water contents with respect to depth and the conditions at the boundaries of the flow domain need to be specified.

Copyright© Markus Tuller and Dani Or2002-2004