## Biomechanics Of Surface Runoff And Soil Water Percolation

ABSTRACT

In this study, the complex interaction of surface runoff with the biomechanics of soil water transport and heat transfer rate is theoretically investigated using mathematical model that rely on the two phase flows of an incompressible Newtonian fluid (stormwater) within the soil (porous medium) and on the soil surface (runoff).

The flow and heat transfer characteristics within the soil are determined numerically based on Darcy-Brinkman-Forchheimer model for porous medium coupled with appropriate energy equation while analytical approach is employed to tackle the model for interacting surface runoff stormwater. The effects of various embedded biophysical parameters on the temperature distribution and water transport in soils and across the surface runoff together with soil-runoff interface skin friction and Nusselt number are display graphically and discussed quantitatively. It is found that an increase in surface runoff over tightly packed soil lessens stormwater percolation rate but enhances both soil erosion and heat transfer rate.

TABLE OF CONTENTS

Declaration of Authorship iii

Abstract vii

Acknowledgements ix

1 Background of Study 1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Definition of Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1 Fluid [26] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

Fluid Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2.2 Heat Transfer [19] . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Conduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2.3 Groundwater . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Aquifers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2.4 infiltration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2.5 Percolation Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2.6 Soil Erosion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Causes of Soil Erosion . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2.7 Stormwater . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.2.8 Flow in Porous Media . . . . . . . . . . . . . . . . . . . . . . . . 10

1.3 Aims and Objectives of Study . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4 Significance of Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.5 Research Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5.1 Continuity Equation . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5.2 Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . 13

1.5.3 Energy Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.5.4 Shooting method [14] . . . . . . . . . . . . . . . . . . . . . . . . 17

1.5.5 Runge-Kutta Method [12] . . . . . . . . . . . . . . . . . . . . . 18

1.5.6 Runge-Kutta-Fehlberg Method (RKF45) [16] . . . . . . . . . . . 19

1.6 Study Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2 Literature Review 23

3 Model Problem and Solution Procedure 25

3.1 Model Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2 Solution Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

xii

4 Discussion and Graphical Results 29

4.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.1.1 Runoff and Soil Percolation Velocity Profiles . . . . . . . . . . . 29

4.1.2 Surface Runoff and Soil Percolation Temperature Profiles . . . . 29

4.1.3 Soil Surface Erosion (Skin Friction) and Heat Transfer Rate (Nusselt

Number) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.2 Graphical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5 Conclusion 37

CHAPTER ONE

Background of Study

1.1 Introduction

Surface runoff and soil water percolation are closely associated with rainfall and melting of snow, or glaciers. Soil inability to absorb excess stormwater and meltwater due to heavy rainfall, high melt rate of snow and glacier, soil saturation, impervious resulting from surface sealing or pavement, etc., do lead to surface runoff [10]. Surface runoff is the major cause of soil erosion and surface water pollution. In urban areas, runoff is the main cause of flooding which may damage properties and infrastructures including loss of life [23]. In order to alleviate the unpleasant effects of surface runoff, several proactive measures are needed to boost soil absorption of stormwater and meltwater. These measures may include minimizing impervious surfaces in urban areas, adopting soil erosion and flood control programs, etc.

Moreover, percolation describes the downward flow rate of the stormwater or meltwater within the soil [6]. Water percolation in the soil contributes to the formation of groundwater aquifers which serves as a freshwater storage that can be utilized during droughts when surface water supplies are reduced. Generally, soil is regarded as a porous media; the soil loose sediments like sand and gravel are porous and permeable.

It can hold water and allows water to flow through [5]. While the amount of porosity in a soil depends on its mineral content and structure, the rate of water percolation depends on soil permeability (i.e. the size of the soil pore spaces and how the pores are connected). For instance, sandy soils have large well connected pores and higher permeability than the clay soils [17]. The use of mathematical models to tackle the menace of surface runoff and enhance the soil water percolation for the formation of groundwater aquifers has attracted the attention of several scientists and researchers [7, 22, 4, 21, 20, 15]. For soil with weak permeability, the relationship between the flow rates and the pressure gradient would be practically linear based on the Brinkman form of Darcy law, while this relationship may be nonlinear for soil with strong permeability (Darcy-Forchheimer law) [25, 2, 18]. Bristow and Horton [3] theoretically investigated the influence of surface mulch soil water flow and heat transfer. The effects of temperature gradient on the soil water flow were studied by Gurr et al. [9]. Numerical results on soil water flow and heat transfer rate together with soil-atmosphere interaction was reported by Fetzer et al. [8]. In all the above studies, it is observed that mathematical model of soil-runoff interface at the continuum scale where water and energy fluxes are highly dynamic are often

magnified. This may lead to inaccuracy in the result obtained.

In this present study, the biomechanics of surface runoff and its interaction with soil water percolation is numerically examined. The Darcy-Brinkman-Forchheimer nonlinear model for porous medium coupled with appropriate energy equation is employed in order to analysis the soil percolation rate while the model representing.

2 Chapter 1. Background of Study

Interacting surface runoff is based on modified Blasius flow with heat transfer characteristics.

The groundwater aquifer servers as the lower boundary of the porous medium domain while the soil surface represents the upper boundary to the porous medium domain and is dramatically influenced by changes in velocity and temperature gradients of both runoff and stormwater percolation. In chapter 3, the model and its mathematical equations are obtained, analysed and solved.In chapter 4 Pertinent results are graphically presented and discussed. The thesis provides a mathematical treatment of surface runoff menace and veritable platform to understand the complex interaction between the surface runoff and the soil water percolation rate.

1.2 Definition of Terms

In this section we are going to define the important terminologies.

1.2.1 Fluid [26]

A fluid is a substance that deforms continuously when acted on by a shearing stress of any magnitude.

Fluid can be divided into three part:

1. Liquids such as water,oil, and gasoline

2. Gases such as propane, methen, co2 etc.

3. Plasma such as blood plasma.

Fluid Properties

To study any fluid’s behavior, it is necessary to know and discuss some fluid properties such as Density (r), Specific Weight (g), Specific Gravity (SG), and Specific Volume (n).

1. Density (r):

Density is the mass of a fluid per unit volume, and its SI unit is (kilogram/meter3)

Density, r =

Mass

Volume

kg

m3 (1.1)

In general, density of a fluid decreases with increase in temperature. It increases with increase in pressure. The ideal gas equation is given by:

PV = mRT [whereR ! Universal Gas Constan] ,

P =

m

V

RT

P = rRT

h

since, r =

m

V

i

. (1.2)

Equation (1.2) is used to find the density of any fluid, if the pressure (P) and temperature

(T) are known.

Note: The density of standard liquid (water) is 1000 kg/m3

2. Specific Weight (g):

Specific weigth is the weight of a fluid per unit volume, and its SI unit is N/m3

(Newton per meter cubed)

g =

weight

volume

N

m3 (1.3)

1.2. Definition of Terms 3

Also, specific weight is a function of density as described by the following relationship:

g = rg (1.4)

where: g is the acceleration due to gravity, and g = 9.81 (m/s2)

Specific weight varies from place to place due to the change of acceleration due to

gravity (g).

3. Specific Gravity (SG):

Specific gravity is ratio of the density of a fluid to the density of water at (4C), at

this temperature the density of water is 1000 kg/m3 in SI unit.

SG =

rFluid

rH2O at 4C

(1.5)

We can also use the ratio of specific weight of a fluid to the specific weight of water at 4C, at this temperature the specific weight of water is 9.81 KN/m3 (kilo newton per meter cubed) in SI units.

Specific gravity may also be defined as the ratio of specific weight of the given fluid to the specific weight of standard fluid.

SG =

SpecificWeight of Given Fluid

SpecificWeight of Standard Fluid

=

gFluid

gH2O at 4C

(1.6)

4. Specific Volume (n):

Specific volume is the volume of a fluid per unit mass. It is the reciprocal of density.

n =

1

r

m3

kg

(1.7)

Viscosity

Is the term used to describe the fluidity of a fluid.

Dynamic Viscosity

Is the shear stress required to cause a unit change in the rate of angular deformation

of the fluid.

Dynamic viscosity is a fluid property that relates shearing stress (t) to fluid motion

( du

dy ).

By Newton’s law of viscosity.

t = m

du

dy

Thus m = t

dy

du

(1.8)

t = shear stress or FA

force per unit area N

m2

m = dynamic viscosity N.S

m2

du

dy = velocity gradiant or rate of angular deformation (1S

).

Notice that the dynamic viscosity (m) is the most important factor in the equation which is able to control shear stress (t) and gradient ( du dy ).

The relationship between t and du dy is linear with the slope equal to viscosity (m).

Kinematic Viscosity

Is the relationship between dynamic viscosity (m) and density (r) when the force

dimension cancels.

n =

m

r

(1.9)

4 Chapter 1. Background of Study

FIGURE 1.1: Retionship between t and du dy

Where:

n = Kinematic Viscosity (m2

s )

m = Dynamic Viscosity ( N.s

m2 )

r = Density ( kg

m3 ).

Temperature

It is the property that determines the degree of hotness or coldness or the level of heat intensity of a fluid. Temperature is measured by using temperature scales.There are 3 commonly used temperature scales. They are

1. Celsius (or centigrade) scale

2. Fahrenheit scale

3. Kelvin scale (or absolute temperature scale)

Kelvin scale is widely used in engineering. This is because, this scale is independent of properties of a substance.

Pressure

Pressure of a fluid is the force per unit area of the fluid. In other words, it is the ratio of force on a fluid to the area of the fluid held perpendicular to the direction of the force.

Pressure is denoted by the letter ‘P’. Its unit is N/m2.

1.2.2 Heat Transfer [19]

Heat transfer which is defined as the transmission of energy from one region to another as a result of temperature gradient takes place by the following three modes:

(i) Conduction ;

(ii) Convection ;

(iii) Radiation.

Heat transmission, in majority of real situations, occurs as a result of combinations of these modes of heat transfer. Example : The water in a boiler shell receives its heat from the fire-bed by conducted, convected and radiated heat from the fire to the shell, conducted heat through the shell and conducted and convected heat from the inner shell wall, to the water. Heat always flows in the direction of lower temperature.

1.2. Definition of Terms 5

The above three modes are similar in that a temperature differential must exist and the heat exchange is in the direction of decreasing temperature ; each method, however, has different controlling laws.

FIGURE 1.2: Illustration of conduction, convection and radiation heat transfer Conduction

Conduction is the transfer of heat from one part of a substance to another part of the same substance, or from one substance to another in physical contact with it, without appreciable displacement of molecules forming the substance.

In solids, the heat is conducted by the following two mechanisms:

(i) By lattice vibration (The faster moving molecules or atoms in the hottest part of a body transfer heat by impacts some of their energy to adjacent molecules).

(ii) By transport of free electrons (Free electrons provide an energy flux in the direction of decreasing temperature—For metals, especially good electrical conductors, the electronic mechanism is responsible for the major portion of the heat flux except at low temperature).

In case of gases, the mechanisam of heat conduction is simple. The kinetic energy of a molecule is a function of temperature. These molecules are in a continuous random motion ex-changing energy and momentum. When a molecule from the high temperature region collides with a molecule from the low temperature region, it loses energy by collisions.

In liquids, the mechanism of heat is nearer to that of gases. However, the molecules are more closely spaced and intermolecular forces come into play.

Convection

Convection is the transfer of heat within a fluid by mixing of one portion of the fluid with another.

1 Convection is possible only in a fluid medium and is directly linked with the transport of medium itself.

2 Convection constitutes the macroform of the heat transfer since macroscopic particles of a fluid moving in space cause the heat exchange.

3 The effectiveness of heat transfer by convection depends largely upon the mixing motion of the fluid.

6 Chapter 1. Background of Study

This mode of heat transfer is met with in situations where energy is transferred as heat to a flowing fluid at any surface over which flow occurs. This mode is basically conduction in a very thin fluid layer at the surface and then mixing caused by the flow. The heat flow depends on the properties of fluid and is independent of the properties of the material of the surface. However, the shape of the surface will influence the flow and hence the heat transfer.

Free or natural convection: Free or natural convection occurs where the fluid circulates by virtue of the natural differences in densities of hot and cold fluids; the denser portions of the fluid move downward because of the greater force of gravity, as compared with the force on the less dense.

Forced convection: When the work is done to blow or pump the fluid, it is said to be forced convection.

Radiation

Radiation is the transfer of heat through space or matter by means other than conduction or convection. Radiation heat is thought of as electromagnetic waves or quanta (as convenient) an emanation of the same nature as light and radio waves.

All bodies radiate heat ; so a transfer of heat by radiation occurs because hot body emits more heat than it receives and a cold body receives more heat than it emits. Radiant energy (being electromagnetic radiation) requires no medium for propagation and will pass through a vacuum.

Note: The rapidly oscillating molecules of the hot body produce electromagnetic waves in hypothetical medium called ether. These waves are identical with light waves, radio waves and X-rays, differ from them only in wavelength and travel with an approximate velocity of 3 108 m/s. These waves carry energy with them and transfer it to the relatively slow-moving molecules of the cold body on which they happen to fall. The molecular energy of the later increases and results in a rise of its temperature. Heat travelling by radiation is known as radiant heat.

The properties of radiant heat in general, are similar to those of light. Some of the properties are:

(i) It does not require the presence of a material medium for its transmission.

(ii) Radiant heat can be reflected from the surfaces and obeys the ordinary laws of reflection.

(iii) It travels with velocity of light.

(iv) Like light, it shows interference, diffraction and polarisation etc.

(v) It follows the law of inverse square.

The wavelength of heat radiations is longer than that of light waves, hence they are invisible to the eye.

1.2.3 Groundwater

Groundwater is fresh water in the rock and soil layers beneath Earth’s land surface.

Some of the precipitation (rain, snow, sleet, and hail) that falls on the land soaks into

Earth’s surface and becomes groundwater. Water-bearing rock layers called aquifers

are saturated (soaked) with groundwater that moves, often very slowly, through

small openings and spaces. This groundwater then returns to lakes, streams, and

marshes (wet, low-lying land with grassy plants) on the land surface via springs

and seeps (small springs or pools where groundwater slowly oozes to the surface).

Groundwater makes up more than one-fifth (22%) of Earth’s total fresh water supply,

1.2. Definition of Terms 7

and it plays a number of critical hydro-logical (water-related), geological and biological

roles on the continents. Soil and rock layers in groundwater recharge zones (a entry

point where water enters an aquifer) reduce flooding by absorbing excess runoff

after heavy rains and spring snowmelts. Aquifers store water through dry seasons

and dry weather, and groundwater flow carries water beneath arid (dry) deserts

and semi-arid grasslands. Groundwater discharge replenishes streams, lakes, and

wetlands on the land surface and is especially important in arid regions that receive

limited rainfall. Flowing groundwater interacts with rocks and minerals in aquifers,

and carries dissolved rock-building chemicals and biological nutrients. Vibrant communities

of plants and animals (ecosystems) live in and around groundwater springs

and seeps.

Almost all of the fresh liquid water that is readily available for human use comes

from underground. (The bulk of Earth’s fresh water is frozen in ice in the North

and South Pole regions. Water in streams, rivers, lakes, wetlands, the atmosphere,

and within living organisms makes up only a tiny portion of Earth’s fresh water.)

For thousands of years, humans have used groundwater from springs and shallow

wells to fill drinking water reservoirs, and water livestock and crops. Today, human

water needs far exceed surface water supplies in many regions, and Earth’s rapidlygrowing

human population relies heavily upon groundwater to meet its ever larger

demand for clean, fresh water.

FIGURE 1.3: Groundwater

Aquifers

An aquifer is a body of rock or soil that yields water for human use. Most aquifers

are water-saturated layers of rock or loose sediment. With the exception of a few

aquifers that have water-filled caves within them, aquifers are not underground

lakes or holding tanks, but rather rock “sponges” that hold groundwater in tiny

cracks, cavities, and pores (tiny openings in which a liquid can pass) between mineral

grains (rocks are made of minerals). The total amount of empty pore space in the

rock material, called its porosity, determines the amount of groundwater the aquifer

can hold. Materials like sand and gravel have high porosity, meaning that they can

absorb a high amount of water. Rocks like granite, marble, and limestone have low

porosity, and make poor groundwater reservoirs.

Aquifers must have high permeability in addition to high porosity. Permeability is

the ability of the rock or other material to allow water to pass through it. The pore

space in permeable materials is interconnected throughout the rock or sediment, allowing

groundwater to move freely through it. Some high-porosity materials, like

mud and clay, have very low permeability. They soak up and hold water, but don’t

release it easily to wells or other groundwater discharge points, so they are not good

aquifer materials. Sandstone, limestone, fractured granite, glacial sediment, loose

sand, and gravel are examples of materials that make good aquifers.

8 Chapter 1. Background of Study

FIGURE 1.4: Aquifer

1.2.4 infiltration

Infiltration is the process by which water on the ground surface enters the soil.

Infiltration is governed by two forces, gravity, and capillary action.

While smaller pores offer greater resistance to gravity, very small pores pull water

through capillary action in addition to and even against the force of gravity.

Infiltration rate in soil science is a measure of the rate at which a particular soil is

able to absorb rainfall or irrigation.

It is measured in inches per hour or millimeters per hour.

The rate decreases as the soil becomes saturated.

If the precipitation rate exceeds the infiltration rate, runoff will usually occur unless

there is some physical barrier.

It is related to the saturated hydraulic conductivity of the near-surface soil.

FIGURE 1.5: Infiltration

1.2.5 Percolation Rate

When we sprinkle water on the ground, it is soon absorbed by the soil. This is

because water percolates through the soil. The process in which water passes down

slowly through the soil is called percolation of water. But water does not percolate

at the same rate in all types of soils.

Sandy soil allows maximum percolation of water and clay soil allows minimum

percolation of water. Rainwater percolates through the soil and collects above the

bedrock. This level of groundwater is called water table. Sandy soil is quite loose,

so the percolation rate of water is highest in sandy soil but lowest in the clay soil

because it is very compact.

Paddy (rice crops) is planted in standing water in the fields. Hence, the soil with a

low percolation rate of water would be the most suitable for growing paddy because

it will allow the water to remain in the fields for a much longer time.

1.2. Definition of Terms 9

1.2.6 Soil Erosion

t is a process in which the top fertile layer of soil is lost. Due to soil erosion, the soil

becomes less fertile. The top layer of soil is very light which is easily carried away

by wind and water. The removal of topsoil by the natural forces is known as soil

erosion.

Causes of Soil Erosion

Various agents, like wind, water, deforestation, overgrazing by cattle, etc., cause soil

erosion. The various factors of soil erosion are:

1. Wind

When strong winds blow, the topsoil along with the organic matter is carried away

by the wind. This happens more often when the land is not covered with grass or

plants. Such conditions are very common in desert and semi-desert regions where

strong winds blow very frequently.

2. Water

When it rains in the hilly areas, the soil gets washed away towards the plains. The

running water deposits the mineral-rich soil in the riverbed and over the years this

deposition of soil can change the course of the river. This can lead to floods which

cause the destruction of life and property.

FIGURE 1.6: Water erosion leads to loss of agricultur potential

3. Overgrazing

When cattle are allowed to graze on the same field repeatedly, all the available grass,

including the roots are eaten by them. This makes the topsoil vulnerable to wind and

flowing water, leading to soil erosion.

4. Deforestation

Humans have taken land from the forest to cultivate in order to feed the everincreasing

population and to build houses, industries, etc. Cutting down of trees

on a large scale for these purposes is deforestation. The roots of trees hold the soil

together, thus preventing the soil from getting uprooted. When large areas of the

forest are cleared, the topsoil gets eroded by wind and flowing water.

1.2.7 Stormwater

Stormwater is a term used to describe water that originates during precipitation

events. It may also be used to apply to water that originates with snowmelt or runoff

10 Chapter 1. Background of Study

Water from overwatering that enters the stormwater system. Stormwater that does

not soak into the ground becomes surface runoff, which either flows directly into

surface waterways or is channeled into storm sewers, which eventually discharge to

surface waters and the ocean.

Stormwater is of concern for two main issues: one related to the volume and timing

of runoff water (flood control and water supplies) and the other related to potential

contaminants that the water is carrying, i.e. water pollution. Stormwater outfalls

that discharge directly on a beach often create shallow streams or pools of water that

are an “attractive nuisance.” This water often has high concentrations of fecal indicator

bacteria and may contain human pathogens. Contact with such water, especially

by small children, should be avoided.

FIGURE 1.7: South Carolina stormdrain outfall and warning sign

1.2.8 Flow in Porous Media

A porous medium can be defined as a material consisting of solid matrix with an

interconnected void. The interconnected pores are very important because they are

the ones that affect the flow. The definition of the porosity of the porous medium

can be given as the ratio of pore volume to the total volume of a given sample of

material. Darcy’s law [18] expressed as

rP =

m

K

Q (1.10)

is used for flow in porous media, where K is the porous medium permeability, is

the fluid dynamic viscosity, P is the fluid pressure and Q is the volumetric flow

rate. It works with variables averaged over several pore widths. Darcy’s law may

be extended to include transitional flow between boundaries (i.e. Darcy -Brinkman

model),

mr2q rP

mf

K

q = 0, (1.11)

where mf is known as the effective viscosity and q is the velocity vector. In order to

account for the nonlinear behaviour of the pressure difference, the inertial term due

to Forchheimer [18] may be add (i.e Darcy -Brinkman-Forchheimer model)

mr2q rP

mf

K

q

cr

p

K

q2 = 0 (1.12)

1.3. Aims and Objectives of Study 11

where r is the fluid density and c is the porous medium inertial parameter. In recent

years, the investigation of flow of fluids through porous media has become an

important topic due to the recovery of crude oil from the pores of reservoir rocks.

Also the flow through porous media is of interest in chemical engineering (absorption,

filtration), petroleum engineering, hydrology, soil physics, bio-physics and geophysics.

1.3 Aims and Objectives of Study

The main objectives of the study in this thesis are as follows:

To derive a nonlinear mathematical model for porous medium coupled with appropriate energy equation.

To solve the two points boundary value problem of soil water percolation numerically using Runge-Kutta-Fehlberg integration scheme coupled with shooting method.

To illustrate the impact of embedded biophysical parameters on stormwater velocity profiles both within the soil in the region 0 h 1 and the runoff in the region h > 1.

To show the effects of various biophysical parameters on the stormwater temperature profiles both within the soil in the region 0 h 1 and the runoff in the region h > 1.

To demonstrate the effects of various biophysical parameters on the coefficient of skin friction which invariably lead to soil erosion and the heat transfer rate at the soil surface (h = 1) due to interaction between the runoff and the soil water percolation.

1.4 Significance of Study

Surface runoff is a water from rain, snowmelt, or other sources, that flows over the land surface, and is a major component of the water cycle. When runoff flows along the ground, it can pick up soil contaminants such as petroleum, pesticides, or fertilizers that become discharge or overland flow. Urbanization increases surface runoff, by creating more impervious surfaces such as pavement and buildings do not allow percolation of the water down through the soil to the aquifer. Increased runoff reduces groundwater recharge, thus lowering the water table and making droughts worse, especially for farmers and others who depend on water wells. Runoff is an economic threat, as well as an environmental one. Agribusiness loses millions of dollars to runoff every year. In the process of erosion, runoff can carry away the fertile layer of topsoil. Farmers rely on topsoil to grow crops. Tons of topsoil are lost to runoff every year. People can limit runoff pollution in many ways. Farmers and gardeners can reduce the amount of fertilizer they use. Urban areas can reduce the number of impervious surfaces. Soil acts as a natural sponge, filtering and absorbing many harmful chemicals. Communities can plant native vegetation. Shrubs and other plants prevent erosion and runoff from going into waterways.

12 Chapter 1. Background of Study

1.5 Research Methodology

1.5.1 Continuity Equation

The system is a fixed quantity of mass, denoted by m.Thus the mass of the system is

conserved and does not change.

msyst = const

or

dm

dt

= 0. (1.13)

The total mass of the fluid is:

dm = r dV (1.14)

equation 1.14 implies

m =

Z

dm

=

Z

V

r dV (1.15)

By substituting equation 1.15 into equation 1.13 we obtian:

d

dt

Z

V

r dV = 0 (1.16)

By using Reynold transport theorem:

d

dt

Z

V

r dV =

Z

V

¶r

¶t

dV +

Z

S

r (q. ˆn)dS = 0 (1.17)

Now by using Gauss divergence theorem onto

R

S

r(q. ˆn)dS we have:

Z

V

¶r

¶t

dV +

Z

V

r.(rq) dV = 0 (1.18)

)

Z

V

¶r

¶t

+ r.(rq)

dV = 0 (1.19)

)

¶r

¶t

+ r.(rq) = 0 (1.20)

Equation 1.20 is called the continuity equation.

Where:

r is the mass density

q is the velocity vector.

when the density (r) is constant equation 1.20 becomes r.q = 0 and this is the case

for incompressible flow.

1.5. Research Methodology 13

1.5.2 Momentum Equation

Momentum (p) is the product of mass (m) and velocity (q):

p = mq (1.21)

Note that q(u, v,w) is velocity vector and recall that:

r =

m

V

) m = rV (1.22)

Where: r is the density, m is the mass, and V is the volume.

By substituting equation 1.22 into equation 1.21 we obtain:

p = mq = (rq)V = qrV

) dp = rq dV (1.23)

By taking integral of equation 1.23 in both sides we obtain:

Z

dp =

Z

rq dV

) p =

Z

V

rq dV (1.24)

By differentaiting equation 1.24 with respect to time (t), we have:

dp

dt

=

d

dt

Z

V

rq dV = F

Note that: Force (F) is the rate of change of momentum with respect to time (t).

let us assume an incompressible fluid (i.e. r is constant):

d

dt

Z

V

rq dV = Fs + Fb (1.25)

Where: Fs is the surface force, and Fb is the body force.

Fs =

Z

S

t. ˆn dS

Fb =

Z

V

r fm dV (1.26)

Then substituting equation 1.26 into equation 1.25 we have:

d

dt

Z

V

rq dV =

Z

S

t. ˆn dS +

Z

V

r fm dV (1.27)

By using Reynold transport theorem equation 1.27 becomes:

Z

V

¶

¶t (pq) dV +

Z

S

rq(q.nˆ) dS =

Z

S

t. ˆn dS +

Z

V

r fm dV (1.28)

14 Chapter 1. Background of Study

Using Gauss divergence theorem equation 1.28 becomes:

Z

V

¶

¶t (pq) dV +

Z

V

r.(rq)q dV =

Z

V

r.t dV +

Z

V

r fm dV

)

Z

V

r

¶q

¶t

+ r.(rq)q r.t r fm

dV = 0

) r

¶q

¶t

+ r.(rq)q = r.t + r fm (1.29)

But r.(rq)q = (rq.r)q + rq(r.q), for incompressible fluid r.q = 0, therefore,

r.(rq)q = (rq.r)q (1.30)

By substituting equation 1.30 into equation 1.29 we obtain:

r

¶q

¶t

+ (rq.r)q = r.t + r fm (1.31)

t =

0

@

p + txx txy txz

tyx p + tyy tyz

tzx tzy p + tzz

1

A = p

0

@

1 0 0

0 1 0

0 0 1

1

A +

0

@

txx txy txz

tyx tyy tyz

tzx tzy tzz

1

A

Where:

txx = 2m¶u

¶x , tyy = 2m¶v

¶y , tzz = 2m¶w

¶z

txy = tyx = m

¶u

¶y + ¶v

¶x

tyz = tzy = m

¶w

¶y + ¶v

¶z

txz = tzx = m

¶w

¶x + ¶u

¶z

r.t =

¶p

¶x

+ mr2u

i +

¶p

¶y

+ mr2v

j +

¶p

¶z

+ mr2w

k

= rp + mr2q (1.32)

By substituting equation 1.32 into equation 1.31 we get:

r

¶q

¶t

+ (rq.r)q = rp + mr2q + r fm (1.33)

By dividing both sides of equation 1.33 by r we have:

¶q

¶t

+ (q.r)q =

1

r

rp +

m

r

r2q + fm (1.34)

Let n = m

r and substitute it into equation 1.34 to get:

¶q

¶t

+ (q.r)q =

1

r

rp + nr2q + fm (1.35)

Equation 1.35 is called Navier-Stokes equation for an incompressible viscous fluid.

Where:

q is the velocity vector of the flow.

p is the pressure.

1.5. Research Methodology 15

fm is the body force.

n is the kinematic viscosity.

1.5.3 Energy Equation

From the first law of thermodynamics DE = W + Q where E is the total internal

energy of the system, W is the work done, and Q is the heat transfer.

Note that: Q can be positive or negative.

dE

dt

=

dQ

dt

+

dW

dt

(1.36)

We want to use equation 1.36 to get the energy equation for fluid.

Total energy of the system = kinetic energy + potential energy.

dE = redV +

r

2

q2dV where e = cvt = cpT

The total energy is E =

R

re dV +

R rq2

2 dV therefore,

dE

dt

=

d

dt

Z

re dV +

Z

rq2

2

dV

dE

dt

=

d

dt

Z

V

re +

rq2

2

dV (1.37)

dQ

dt

=

Z

S

(H. ˆn) dS (1.38)

Where H is the heat transfer rate.

Fourer Heat Conduction Law gives us an expression for heat transfer which is given

as H = krT where k is thermal conductivity, also from equation 1.36 dW

dt = work

done by surface force (qFs) + work done by body force (qFb)

dW

dt

=

Z

S

q(t. ˆn) dS +

Z

V

r(q fm) dV (1.39)

Let’s substitute equations (1.37, 1.38, and 1.39) into equation 1.36 to obtain:

d

dt

Z

V

re +

rq2

2

dV =

Z

S

(krT. ˆn) dS +

Z

S

q(t. ˆn) dS +

Z

V

r(q fm) dV

We applied Reynolds Transport theorem:

Z

V

¶

¶t

re +

rq2

2

dV +

Z

S

q

re +

rq2

2

. ˆn

dS =

Z

S

krT. ˆn dS +

Z

S

q(t. ˆn) dS +

Z

V

r(q fm) dV

)

Z

V

¶re

¶t

dV +

Z

V

¶

¶t

rq2

2

dV +

Z

S

re(q. ˆn) dS +

Z

S

rq2

2

(q. ˆn) dS =

Z

S

krT. ˆn dS +

Z

S

q(t. ˆn) dS +

Z

V

r(q fm) dV

16 Chapter 1. Background of Study

Therefore, we have

Z

V

¶re

¶t

dV +

Z

V

rq ¶q

¶t

dV +

Z

S

re(q. ˆn) dS +

Z

S

rq2

2

(q. ˆn) dS =

Z

S

krT. ˆn dS +

Z

S

q(t. ˆn) dS +

Z

V

r(q fm) dV

Now let’s apply Gauss Divergence Theorem:

Z

V

r

¶e

¶t

dV +

Z

V

rq ¶q

¶t

dV +

Z

V

r.(req) dV +

Z

V

r.(

rq2

2

q) dV =

Z

V

r.(krT) dV +

Z

V

q(r.t) dV

+

Z

V

r(q fm) dV (1.40)

Note r. rq2

2 q = q.r( rq2

2 ) + rq2

2 r.q (r.q = 0 for incompressible)

=) q.r( rq2

2 ) = q. ¶

¶q ( rq2

2 ) ¶q

¶x = q.(rq)rq

In the same way r.qt = q(r.t) + (t.r)q

Now equation 1.40 can be written as:

Z

V

r

¶e

¶t

dV +

Z

V

rq ¶q

¶t

dV +

Z

V

r.(req) dV +

Z

V

q(rq)rq dV =

Z

V

r.(krT) dV +

Z

V

q(r.t) dV

+

Z

V

(t.r)q dV +

Z

V

rq fm dV

(1.41)

Note that: By using Navier-Stokes Equation

R

V

rq

¶q

¶t + q.rq

dV =

R

V

q (r.t + r fm) dV

then equation 1.41 becomes: R

V

r ¶e

¶t + r(r.eq)

dV =

R

V

(r.krT + (t.r)q) dV

This implies that r

¶e

¶t + r.eq

= r.krT + (t.r)q

Note: r.eq = e(r.q) + (q.r)e

Therefore, r

¶e

¶t + e(r.q) + (q.r)e = r.krT + (t.r)q

, and (r.q = 0)

But e = cpT, therefore:

rcp

¶T

¶t

+ q.rT

= r.krT + (t.r)q

If k is constant then we have:

rcp

¶T

¶t

+ q.rT

= kr2T + (t.r)q

Note: (t.q)q = F (viscous dissipation function)

rcp

¶T

¶t

+ q.rT

= kr2T + F

¶T

¶t

+ q.rT =

k

rcp

r2T +

1

rcp

F (1.42)

1.5. Research Methodology 17

Equation 1.42 is called Energy Equation.

Where:

q.rT is the heat convection term.

k

rcp

r2T is the heat conduction term.

1

rcpF is the dissipation term.

Note that: In cartesian coordinates

F = m

”

2

¶u

¶x

2

+ 2

¶v

¶y

2

+ 2

¶w

¶z

2

+

¶v

¶x

+

¶u

¶y

2

+

¶w

¶y

+

¶v

¶z

2

+

¶u

¶z

+

¶w

¶x

2

#

1.5.4 Shooting method [14]

Consider a 2nd order ordinary differential equation with two boundary conditions

y00 = f (x, y, y0), a < x < b

y(a) = a

y(b) = b,

where a, b, a, b are given constants, y is the unknown function of x, f is a given

function that specifies the differential equation. This is a two-point boundary value

problem. An initial value problem (IVP) would require that the two conditions be

given at the same value of x. For example, y(a) = a and y0(a) = g. Because the two

separate boundary conditions, the above two-point boundary value problem (BVP)

is more difficult to solve.

The basic idea of “shooting method” is to replace the above BVP by an IVP. But of

course, we do not know the derivative of y at x = a. But we can guess and then

further improve the guess iteratively. More precisely, we treat y0(a) as the unknown,

and use secant method or Newton’s method (or other methods for solving nonlinear

equations) to determine y0(a).

We introduce a function u, which is a function of x, but it also depends on a parameter

t. Namely, u = u(x; t). We use u0 and u00 to denote the partial derivative of u,

with respect to x. We want u to be exactly y, if t is properly chosen. But u is defined

for any t, by

u00 = f (x, u, u0)

u(a; t) = a

u0(a; t) = t.

If you choose some t, you can then solve the above IVP of u. In general u is not the

same as y, since u0(a) = t 6= y0(a). But if t is y0(a), then u is y. Since we do not know

y0(a), we determine it from the boundary condition at x = b. Namely, we solve t

from:

f(t) = u(b; t) b = 0.

If a solution t is found such that f(t) = 0, that means u(b; t) = b. Therefore, u

satisfies the same two boundary conditions at x = a and x = b, as y. In other words,

u = y. Thus, the solution t of f(t) = 0 must be t = y0(a).

If we can solve the IVP of u (for arbitrary t) analytically, we can write down a formula

for f(t) = u(b; t) b. Of course, this is not possible in general. However, without

an analytic formula, we can still solve f(t) = 0 numerically. For any t, a numerical

18 Chapter 1. Background of Study

method for IVP of u can be used to find an approximate value of u(b; t) (thus f(t)).

The simplest method is to use the secant method.

tj+1 = tj

tj tj1

f(tj) f(tj1)

f(tj) , j = 1, 2, 3, . . .

For that purpose, we need two initial guesses: t0 and t1 . We can also use Newton’s

method:

tj+1 = tj

f(tj)

f0(tj)

, j = 0, 1, 2, . . .

We need a method to calculate the derivative f(t). Since f(t) = u(b; t) b, we have

f0(t) =

¶u

¶t

(b; t) 0 =

¶u

¶t

(b; t).

If we define v(x; t) = ¶u

¶t , we have the following IVP for v:

v00 = fu(x, u, u0)v + fu0 (x, u, u0)v0

v(a; t) = 0

v0(a; t) = 1.

Here v0 and v00 are the first and 2nd order partial derivatives of v, with respect to x.

The above set of equations are obtained from taking partial derivative with respect to x for the system for u. The chain rule is used to obtain the differential equation of v. Now, we have f0(t) = v(b; t).

1.5.5 Runge-Kutta Method [12]

The Runge-Kutta method is the most widely used method of solving differential equations with numerical methods. It differs from the Taylor series method in that we use values of the first derivative of f (x, y) at several points instead of the values of successive derivatives at a single point.

For a Runge-Kutta method of order 2, the following formulas are applicable.

k1 = h f (xn, yn)

k2 = h f (xn + h, yn + h)

yn+1 = yn +

1

2

(k1 + k2) (1.43)

Equation 1.43 is for Runge-Kutta Method of order 2.

When higher accuracy is desired, we can use order 3 or order 4. The applicable

formulas are as follows:

k1 = h f (xn, yn)

k2 = h f (xn +

h

2

, yn +

k1

2

)

k3 = h f (xn + h, yn + 2k2 k1)

yn+1 = yn +

1

6

(k1 + 4k2 + k3) (1.44)

1.5. Research Methodology 19

Equation 1.44 is for Runge-Kutta Method of order 3

k1 = h f (xn, yn)

k2 = h f (xn +

h

2

, yn +

k1

2

)

k3 = h f (xn +

h

2

, yn +

k2

2

)

k4 = h f (xn + h, yn + k3)

yn+1 = yn +

1

6

(k1 + 2K2 + 2k3 + k4) (1.45)

Equation 1.45 is for Runge-Kutta Method of order 4

1.5.6 Runge-Kutta-Fehlberg Method (RKF45) [16]

One way to guarantee accuracy in the solution of an I.V.P. is to solve the problem twice using step sizes h and h/2 and compare answers at the mesh points corresponding to the larger step size. But this requires a significant amount of computation for the smaller step size and must be repeated if it is determined that the agreement is not good enough.

The Runge-Kutta-Fehlberg method (denoted RKF45) is one way to try to resolve this problem. It has a procedure to determine if the proper step size h is being used. At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step size is reduced. If the answers agree to more significant digits than required, the step size is increased.

Each step requires the use of the following six values:

k1 = h f (xn, yn)

k2 = h f (xn +

1

4

h, yn +

1

4

k1),

k3 = h f

xn +

3

8

h, yn +

3

32

k1 +

9

32

k2

,

k4 = h f

xn +

12

13

h, yn +

1932

2197

k1

7200

2197

k2 +

7296

2197

k3

,

k5 = h f

xn + h, yn +

439

216

k1 8k2 +

3680

513

k3

845

4104

k4

,

k6 = h f

xn +

1

2

h, yn

8

27

k1 + 2k2

3544

2565

k3 +

1859

4104

k4

11

40

k5

(1.46)

Then an approximation to the solution of the I.V.P. is made using a Runge-Kutta

method of order 4:

yn+1 = yn +

25

216

k1 +

1408

2565

k3 +

2197

4101

k4

1

5

k5 (1.47)

where the four function values k1 , k3 , K4 , and k5 are used. Notice that k2 is not used

in formula (1.47). A better value for the solution is determined using a Runge-Kutta

method of order 5:

zn+1 = yn +

16

135

k1 +

6656

12825

k3 +

28561

56430

k4

9

50

k5 +

2

55

k6. (1.48)

20 Chapter 1. Background of Study

1.6 Study Limitations

The work in the present was done with the help of mathematical tools and methods in which some assumption have been made for an accuracy of the results. For instance, in order to gain an insight into the complex interaction between the surface runoff and soil water percolation, some assumption were made for the values of biophysical parameters that were utilised for the numerical computations. This study only consider the problem on the two phase flows of an incompressible Newtonian fluid within the soil and on the soil surface. It could be extended to the three phase flows of an incompressible Newtonian.