EROSION Submodel

L.J. Hagen

TABLE OF CONTENTS

INTRODUCTION

Submodel Objectives

submodel Outputs

SIMULATION REGION

Spatial Geometry

Input Winds

Wind Barriers

Input Surface Conditions

EROSION THEORY AND ANALYTIC SOLUTIONS

Problem Overview

Saltation/creep Component

Solution for Saltation/creep Discharge

Suspension Component

Solution for Suspension Discharge

PM-10 Component

Solution for PM-10 Discharge

UPDATE SURFACE CONDITIONS DURING WIND EROSION

Surface Update Intervals

Change in Mobile Surface Soil

Update of Crusted Surface

Update of Aggregated Surface

Update of Surface Rock Volume

Update of Soil Fraction < 2.0 mm Diameter

Update of Soil Fraction < 0.10 mm Diameter

Update Surface Roughness

PARAMETERS FOR ANALYTIC SOLUTIONS

Overview

Aerodynamic Roughness

Friction Velocity

Wind Barrier Parameters

Threshold Friction Velocity

Emission Parameters

Abrasion Parameters

Breakage Parameters

Trapping parameters

REFERENCES

INTRODUCTION

The Wind Erosion Prediction System (WEPS) is a continuous, process-based model that simulates weather, field conditions, and wind erosion on crop lands (Hagen, 1991a). The WEPS is composed of a user-interface, databases, and a simulation code. Many of the processes are simulated on a daily basis, but erosion is simulated on a sub-hourly basis. In general, the surface conditions that control wind erosion are all temporal. Hence, the WEPS simulation code is modular in structure and includes weather simulators and five submodels that simulate surface conditions including crop growth, residue decomposition, soil aggregate/crust status, hydrology, and management. On days when wind speed exceeds the threshold for erosion, the erosion submodel then simulates wind erosion over the course of the day. To facilitate estimating onsite and offsite impacts of erosion, the submodel predicts the various components of the soil movement. The EROSION submodel consists of several subroutines written in FORTRAN.

Submodel Objectives

The EROSION submodel objectives are to use parameters supplied by other submodels that describe the soil surface, flat biomass cover, standing biomass leaf and stem areas, wind barriers and weather to:

a) determine if wind erosion can occur in a simulation region on a given day,

b) simulate the processes of wind erosion, when erosion occurs and

c) update changes in the soil surface during erosion events.

Both the surface descriptions and the simulated erosion processes were defined so that each can be experimentally measured (Hagen 2001, Mirzamostafa et al. 1998).

Submodel Outputs

The submodel outputs include estimates of both total soil loss from the simulation region as well as components of the soil loss moving across each field boundary. These components are saltation/creep (> 0.10 mm diameter) and suspended soil (<0.10 mm diameter). The mass of particles less than 10 microns in aerodynamic diameter (PM-10) in the suspended soil also is predicted. Files of time-varying emissions for use as source terms in dust diffusion and transport models can also be output. To compare simulation results to measured erosion data, input files of measured weather can be used in place of the simulated weather.

SIMULATION REGION

Spatial Geometry

In the EROSION submodel, the erodible simulation region is rectangular in shape and surrounded by a boundary that does not allow moving soil in the saltation/creep transport mode to enter the simulation region from upwind. (Fig. E-1). The sides of the erodible simulation region rectangle are assumed oriented parallel to the X and Y axis of a Cartesian coordinate system. Submodel inputs required to define the simulation region include distances from the origin along two coordinate directions to the two diagonal corners (Sim.X1, Y1 and Sim. X2, Y2 ) that are nearest and farthest from the origin, respectively. The rotation angle from north of the erodible simulation region Y- coordinate is also required. The limits on the rotation angle from north are ± 45 degrees.

To further define the region, an automated grid generator is used to overlay a grid of rectangular cells on the erodible simulation region. The erosion is calculated in individual grid cells using a single pass over all the cells. Using a single pass, requires that the sequence of calculations be selected so that they always precedes along the downwind direction. The direction of the calculation passes is only updated daily, because the current wind generator supplies a single wind direction for each day.

In WEPS1.0 the generated grid has a maximum of size of interior 29 X 29 grid cells with a minimum of 7 m cell length. In narrow strips 3 grid cells are always generated, so the minimum cell length may be less than 7 m in this case. The number of grid cells was restricted to reduce computation time, but can easily be modified in the code. When barriers are present, a maximum of 59X59 grid cells are generated. The increased number of grid cells is necessary to permit adequate resolution of the barrier reductions in the wind speed.

The erodible region is assumed to be uniform at the beginning of erosion, but differential erosion among the cells during an erosion event gives rise to a nonuniform surface. The cells are periodically updated during an event to reflect the changes imposed on the surface conditions by erosion. Surface updating enables the model to simulate both source limited and wind energy limited erosion events.

Figure E-1. Simulation region geometry in Cartesian coordinate system. User inputs include opposite corners of the rectangular, erodible, simulation region, the rotation angle from north, and end points of barriers as well as other barrier properties.

Input Winds

Simulated, hourly, wind speeds (m/s) from a weather station representative of the simulation region are normally input to the submodel. The wind speeds generated from the WEPS wind database have been adjusted to a 10 m height (van Donk et al., 2005). The aerodynamic roughness at the weather stations is assumed to be 25 mm. A simulated daily wind direction is also input to the submodel. The EROSION submodel is not called unless maximum daily wind speed at 10 m height reaches 8 m/s. Then, the maximum daily wind speed is used to determine if erosion can occur in the simulation region. If snow depth exceeds 20 mm, no erosion is simulated. Measured wind data at the simulation field or from other heights also may be used as input to the submodel by changing some internal flags.

Wind Barriers

Up to five wind barriers can be placed arbitrarily on erodible areas or boundaries of the simulation region (Fig. E-1). The submodel inputs required to define each wind barrier include the distances to the barrier ends along the two coordinate directions (Bar. X1, Y1 and Bar. X2, Y2 ). Other inputs needed to define the barrier include the height, optical porosity, and total width of the barrier. (The user interface in WEPS 1.0 currently restricts the maximum number of barriers to four and allows them to be placed only on the simulation region boundaries.)

Input Surface Conditions

The surface conditions considered are combinations of the following:

a)	Surface roughness - random and/or oriented measured below the standing biomass canopy but on top of flat residues (Fig. E-2);

b)	Surface covers - flat, random, biomass cover; crust with loose, erodible soil on crust; aggregated soil; and rock cover (> 2.0 mm dia.) (Fig. E-3);

c)	Surface soil moisture, crust/consolidated zone parameters (dry stability and thickness), aggregate parameters (dry stability and size distribution); and

d)	Standing biomass parameters (leaf area index, stalk silhouette area index, height, row spacing, position on ridge, and orientation) (Fig. E-4).

Figure E-2. Illustration of random roughness shelter angles (SA) and oriented roughness spacing (SX-RG) and height (SZ-RG) used in EROSION submodel.

The random roughness parameter is the standard deviation (mm) of the surface roughness below the biomass canopy and includes the roughness contributed by flat biomass.

The oriented roughness parameters include the ridge height, ridge spacing, top width of beds, furrow dike spacing and furrow dike height (mm). The latter are assumed to be the same height as the ridges. Based on the oriented and random roughness parameters, the horizontal fraction of surface area sheltered from saltation impacts is also calculated in the submodel (Potter, Zobeck, and Hagen, 1990). The saltation impact angle was assumed to be 12 degrees above horizontal.

Soil scientists in the U.S. generally report the sum of the soil mass fractions less than 2 mm diameter as 1.0 and then report volume of rocks as a separate value in their databases. As inputs to EROSION, we have followed this precedent and let the numerical values of the surface fractions of crusted and aggregated soil sum to 1.0.

However, in applying surface cover in the EROSION, we assigned rock > 2.0 mm , aggregated soil, and crusted soil as the first level of cover, and forced their effective fractions to sum to 1.0 (Fig. E-3). The second level of cover was assigned as fraction of the crust covered with loose soil; the maximum for this fraction is 1.0. The third level was assigned as the biomass flat fraction of cover, which is assumed to have random distribution over the entire surface and has a maximum value of 1.0.

In addition to the cover fractions illustrated in Fig. E-3, additional input parameters are used to describe the surface soil components. These include loose mass on the crust (kg m-2), crust dry stability [ln(J/kg)], and crust/consolidated zone thickness (mm). Aggregates are described by their dry aggregate stability [ln(J/kg)] (Hagen et al., 1992) and their size distribution using a four parameter modified log-normal distribution function (Wagner and Ding, 1994).

Figure E-3. Diagram of the components of surface fraction (S.F.) cover at each of 3 levels in WEPS.

Standing biomass input parameters include the average height, leaf area index, and stem

area index (Fig. E-4). The latter is defined as the horizontal silhouette area of the stems

per unit ground area. The row spacing, orientation and location (furrow or ridge top) are also input for growing biomass.

Figure E-4. Diagram illustrating above-canopy friction velocity which is reduced by

drag of biomass stems and leaves to the below-canopy friction velocity. The latter is

used to drive the soil erosion processes.

Field surfaces may vary both temporally and spatially. By partitioning complex areas into a series of small, uniform areas and periodically updating the surface conditions, one may thereby encompass both the spatial and temporal variations in fields. However, to simplify user inputs, the initial condition of the simulation region is assumed to have uniform soil and management. To meet stringent run-time requirements, the quasi-steady state solutions to mass conservation equations outlined in this section are used to predict wind erosion.

Soil transport during wind erosion occurs in three modes (Chepil and Woodruff, 1963): creep-size aggregates (0.84 - 2.0 mm diameter) roll along the surface, saltation-size aggregates (0.10 - 0.84 mm diameter) hop over the surface, and suspension-size aggregates ( < 0.10 mm diameter) move above the surface in the turbulent flow. Obviously, as wind speeds, turbulence, or sediment loads change, the diameter of aggregates moving in the various modes also may change slightly (Pye, 1987).

In WEPS, we have assumed that the combined saltation/creep mode of transport has a distinct transport capacity for each surface, based on the surface roughness and wind speed. This assumption generally is supported by both field and wind tunnel measurements of the saltation/creep discharge (Greeley and Iversen, 1985). We also assumed that the suspension component does not reach a transport capacity on most eroding fields. Thus, separate equations were developed for saltation/creep, suspension, and PM10 discharge, because they respond differently to both the wind forces and sediment load (Gillette et al. 1998). Separating these erosion components also is useful, because they have different potential off-site impacts.

Wind erosion occurs over a wide range of surface conditions. To aid in delineating erosion rates among the various surfaces, we identified several individual erosion processes in WEPS (Hagen et al., 1995) (Fig. E-5). These processes include direct entrainment (emission) of loose soil by wind and/or saltation impacts, abrasion of soil from clods/crust by saltation impacts, and breakage of saltation/creep-size aggregates to suspension-size. These processes were selected for individual simulation, because they differ from one another by approximately an order of magnitude in their ability to supply new suspension or saltation/creep-size mass to the airstream in response to a saltation impact (Mirzamostafa et al., 1998). When the saltation/creep discharge exceeds transport capacity over cells in a local area of the surface, trapping of saltation/creep occurs. We also assumed that the coarse fraction of the suspension component begins deposition when moving over cells in the simulation region that are not eroding.

Figure E-5. Diagram illustrating simulated erosion processes on a bare soil in an

on an individual grid cell.

Saltation/creep component

Based on conservation of mass in a control volume (Fig. E-5), a one-dimensional, quasi-
steady state equation for the physical processes involved in saltation/creep is:

where (table 1)

q = horizontal saltation/creep discharge (kg m-1 s-1),

x = downwind distance from nonerodible boundary (m),

Gen = vertical flux from emission of loose aggregates

(kg m-2 s-1),

Gan = vertical flux from abrasion of surface clods and crust

(kg m-2 s-1),

Gssbk = vertical flux of suspension aggregates from

breakage of saltation/creep aggregates (kg m-2 s-1),

Gtp = vertical flux from trapping of saltation/creep

aggregates (kg m-2 s-1),

Each of the vertical fluxes represents either source or sink terms in the control volume and can be estimated by the equations that follow. The net emission source term for loose aggregates is

where

SFssen = mass fraction of suspension-size (<0.10 mm) among loose aggregates (< 2.0 mm diameter),

Cen = coefficient of emission (m-1), and

qen = transport capacity (kg m-1 s-1).

Stout (1990) derived the general form of Eq. 2 and applied it to describe total mass flux at a given height from the surface. However, subsequent research (Hagen, 1991) showed that the abrasion flux from immobile clods and crust was controlled by other factors. Hence, is this report, Eq. 2 applies only to the loose, mobile components of the soil. A typical value for Cen on a loose, bare field is about 0.06 m-1, and values for other conditions have been reported (Hagen et al., 1995).

Many transport capacity equations for saltation/creep have been reported (Greeley and Iversen, 1985). One of the most frequently used was developed by Lettau and Lettau (1978) and can be expressed as

where

Cs = the saltation transport parameter (kg m-4 s2),

with a typical value of about 0.3 or more for surfaces armored with stones,

U* = friction velocity (m s-1), and

U*t = dynamic threshold friction velocity (m s-1).

The suspension-size aggregates are assumed to be mixed intimately with the saltation/creep-size and emitted with them. Although the suspension-size aggregates absorb part of the aerodynamic and impact energy (represented by the emission coefficient) in order to rise from the surface, they do not contribute toward reaching the transport capacity of saltation/creep. Hence, they are subtracted from the total emission of loose aggregates in Eq.2.

The net source term for entrainment of saltation/creep aggregates abraded from immobile clods and crust by impacting saltation/creep is

where

SFssan = mass fraction of suspension-size from abrasion,

Fani = mass fraction saltation impacting clods and crust, and

Cani = coefficient of abrasion (m-1).

An index of two was used in Eq. 4 because, in general, only two targets, exposed clods and crust, must be considered. Other targets, such as residue and rocks, have a Cani near zero. The first term,(1 - Sfssan ), is the fraction of abraded mass that is of saltation/creep-
size. Values of SFssan for some Kansas soils have been measured and ranged from 0.14 to 0.27, depending upon soil texture (Mirzamostafa, 1996). The middle, bracketed term on the right-hand-side of Eq. 4 represents the total soil abraded from clods and crust, as confirmed by wind tunnel experiments (Hagen, 1991). Values for Cani also have been measured for a range of soils and related to their crushing energy (Hagen, et al., 1992). The final term in Eq. 4 is the mass fraction entrained in the air stream. Note that the entrainment rate of this newly created saltation/creep is assumed to be similar to that of loose, saltation/creep-size aggregates already present on the surface, and that the entrainment approaches zero at transport capacity.

A sink for the saltation/creep discharge occurs when these aggregates are broken into suspension-size and carried et al., 1998). This effect is simulated as

where

Cbk = coefficient of breakage (m-1), and

qs = discharge of primary (non-breakable) sand particles (kg m-1 s-1).

The saltation/creep aggregates are more stable than the surface clods and crust, so measured abrasion coefficients average about 9 times more than the breakage coefficients on the same soils (Mirzamostafa, 1996). The wind tunnel experiments also demonstrated that the breakage coefficient remained constant during breakdown of the aggregates to primary particles. The means and variances of these coefficients are related to soil texture. Given q, values for qs can be estimated directly from soil sand content.

Another sink term is the removal of saltation/creep from the air stream by trapping mechanisms (Hagen and Armbrust, 1992). In WEPS, surface trapping and plant interception are simulated as

where

Ct = coefficient of surface trapping (m-1),

Ci = coefficient of plant interception (m-1), and

qcp = transport capacity of the surface, when 40

percent or more is armored (kg m-1 s-1).

When erosive winds cross rough surfaces, such as tillage ridges, that are highly erodible, large amounts of soil are entrained, but a portion of the entrained saltation/creep is often trapped in succeeding downwind furrows. This phenomenon results in a local rearrangement of the surface and reduces net removal of the entrained soil. Our conventionally-defined transport capacity, qen, is based on the threshold velocity where erosion begins. But, when trapping of saltation/creep occurs on rough surfaces, one may hypothesize that qen has been exceeded, and that the true transport capacity of the surface is some value, qcp, that is less than qen. However, qen still appears to be the appropriate limiting value to drive the emission process, because more soil is emitted than can be transported from the local area.

In WEPS, the first term on the right-hand-side of Eq. 6 simulates trapping of saltation/creep by surface roughness. The true transport capacity of the surface, qcp, is based on the threshold friction velocity needed to remove saltation/creep from the furrows. It is calculated using Eq. 3 for a given roughness at the level of clod and crust cover of the surface but with a minimum set at 40 percent of the surface armored. Under this condition, wind tunnel observations show that loose material is removed, but minimal local trapping of saltation/creep occurs.

The second term of Eq. 6 represents interception of saltation/creep by standing plant stalks or other near-surface plant parts. This term arises, because for a given soil surface friction velocity, more transport occurs without than with stalks. This term also is used to assign a higher transport capacity for wind direction parallel to crop rows than for wind direction perpendicular to rows. For saltation normal to the row direction, interception can reduce transport capacity 5 to 10 percent. Comparisons to measured data have been reported previously (Hagen and Armbrust, 1994).

Solution for saltation/creep discharge

When the source and sink terms are collected on the variable q, Eq.1 for saltation/creep can be written in the form

where

and

Integrating eq. 7 along the wind direction, from x1 to x2 and q1 to q2, gives the solution

where

and

Suspension component

Based on conservation of mass in a control volume that extends to the top of the dust cloud, a one-dimensional, quasi-steady state equation for the physical processes generating the suspension component is

or

where

qss = horizontal suspension component discharge (kg m-1 s-1),

Gssen = vertical emission flux of loose, suspension-size aggregates (kg m-2 s-1),

Gssan = vertical flux of suspension-size aggregates created

by abrasion of clods and crust (kg m-2 s-1),

Gssbk = vertical flux of suspension-size aggregates created by breakage of	saltation/creep-size aggregates (kg m-2 s-1),

Cssi = fraction of suspension-size aggregate flux intercepted by standing biomass, and

Gssdp = vertical flux (deposition) of suspension-size aggregates above a non- eroding surface (kg m-2 s-1).

The source and sink terms for the suspension component are simulated by the equations that follow.

For direct emission of loose, suspension-size material by 'splash' impacts and aerodynamic forces

where

Cm = a coefficient of mixing, value about(0.0001 SFssen) (m-1).

Below transport capacity, the driving force causing the emission flux of suspension-size soil is assumed to be similar to that in Eq. E-2 causing the saltation/creep emission flux. This assumption is supported by wind tunnel measurements that show a mixture of suspension-size aggregates and a mixture of saltation-size have about the same threshold velocities (Chepil, 1951).

However, two additional assumptions are inherent in Eq. E-16. The first is that the loose components of saltation/creep and suspension-size aggregates occur as a uniform mixture in the field. As a consequence, during simple net emission, the suspension fraction emitted with the saltation/creep remains the same as it was in the soil. Hence, the suspension fraction can be estimated as

where

SFss = soil mass fraction of loose, suspension-size less than about 0.1 mm, and

SFer = soil mass fraction of loose, erodible-size, less than about 2.0 mm.

The second assumption in Eq. E-16 is that an additional small amount of suspension-size aggregates that are disturbed by the saltation impacts also are entrained, because transport capacity for the suspension component generally is not limiting. The result of this process is gradual depletion of the loose, suspension-size aggregates at the surface. However, when net emission of suspension-size exceeds net emission of saltation/creep-size aggregates, the latter soon dominate the surface area and absorb the impacts, so the process tends to be self-limiting.

For suspension flux created by abrasion of clods and crust,

Additional discussion and measurements of this source term were reported by Mirzamostafa et al.,1998.

For the source of suspension flux created by breakage of saltation/creep aggregates, the term is the same as the sink term in the saltation/creep equation and simulated as

Breakage from impact on immovable targets is assumed to come only from the impacting saltation/creep alone. Breakage coefficients for saltation-size aggregates have been measured in the wind tunnel (Mizamostafa et al.,1998). But the breakage component from impacts on other saltation/creep is assumed to come from both the impacting and target aggregates. Breakage from impact on a mobile target is less likely than breakage from impact on immobile targets. However, these assumptions need further experimental verification.

Finally, the sink term for trapping of suspension flux occurs when the suspension discharge passes over grid cells without active saltation to maintain the suspension flux from the surface. Typically, this implies the presence of a vegetated, water, or rough armored surface. The largest suspension particles, 0.05 to 0.10 mm, comprise roughly half the mass of the suspension discharge (Chepil, 1957; Zobeck and Fryrear, 1986). Through diffusion and settling, they move rapidly toward noneroding surfaces in the simulation region, which serve as sinks. The process is simulated as

where

qsso = maximum value of qss entering deposition region (kg m-1 s-1), and

Cdp = coefficient of deposition of suspension-size (m-1).

The maximum value of Cdp is about 0.02, but less for smooth surfaces or large upwind areas that produce high dust clouds thus, moving a large portion of the soil away from the deposition surface.

Solution for suspension discharge

When the source terms are collected, eq. E-14 can be written in the form:

where

and

Substituting the general solution of eq. E-7, q(x), into eq. E- 21 and integrating along the wind direction from x1 to x2 and qss1 to qss2 give the following equation for suspension discharge

In regions of deposition of suspension component, integration of eq.E- 20 from location x1, with discharge qss1, gives the following for suspension discharge, qss2, at downwind location x2

PM-10 component

Simulation equations for particles with aerodynamic diameters less than 10 m (PM-10 component) of the suspended soil also were developed by subdividing the suspension component emission sources using three parameters to give. (Hagen et al., 1996):

where

SF10en	= soil fraction of PM-10 in suspension-size surface soil,

SF10an = soil fraction of PM-10 in suspension-size aggregates created during abrasion of clods and crust, and

SF10bk	= soil fraction of PM-10 in suspension-size aggregates broken from

saltation and creep-size aggregates.

Solution for PM-10 discharge

When the source terms are collected, eq. E-26 can be written in the form:

where

and

Substituting the general solution of eq. E-7, q(x), into eq. E- 27 and integrating along the wind direction from x1 to x2 and q101 to q102 gives the following equation for PM-10 discharge.

Table 1. List of symbols, definitions, and units.

_________________________________________________________________

Symbols Definitions Units

_________________________________________________________________

Cani	coefficient of abrasion	m-1

Cbk	coefficient of breakage of saltation/creep	m-1

Cdp	coefficient of deposition of suspension-size	m-1

Cen	coefficient of emission	m-1

Ci	coefficient of plant interception	m-1

Cm	coefficient of mixing at soil surface	m-1

Cs	saltation transport parameter	kg m-4 s2

Ct	coefficient of surface trapping	m-1

Fani	mass fraction of saltation impacting clods

and crust

Gan	vertical flux from abrasion of surface clods

and crust	kg m-2 s-1

Gen	vertical flux from emission of loose aggregates	kg m-2 s-1

Gtp	vertical flux from trapping of saltation/creep	kg m-2

aggregates s-1

Gssan	vertical flux of suspension-size aggregates

created by abrasion of clods and crust	kg m-2 s-1

Gssbk	vertical flux of suspension-size aggregates

created by breakage of saltation/creep	kg m-2 s-1

Gssdp	vertical flux (deposition) of suspension-size

aggregates above a non-eroding surface	kg m-2 s-1

Gssen	vertical emission flux of loose,

suspension-size aggregates	kg m-2 s-1

q	horizontal saltation/creep discharge	kg m-1 s-1

qcp	transport capacity of the surface, when 40

percent or more is armored	kg m-1 s-1

qen	transport capacity of saltation/creep	kg m-1 s-1

qs horizontal discharge of primary (non-breakable)

sand particles	kg m-1 s-1

qss	horizontal discharge of suspension component	kg m-1 s-1

qsso	maximum value of qss entering deposition region	kg m-1 s-1

Sfer	soil mass fraction of loose, erodible-size,

less than about 2.0 mm diameter at soil surface

Sfss	soil mass fraction of loose, suspension-size

less than about 0.1 mm diameter at soil surface

SFssan mass fraction of suspension-size (<0.10 mm) of total (<2.0 mm diameter) created by abrasion

Sfssen mass fraction of suspension-size (<0.10 mm)

emitted from loose aggregates (<2.0 mm diameter)

U*	friction velocity	m s-1

U*t	dynamic threshold friction velocity	m s-1

x	downwind distance from nonerodible boundary	m

_________________________________________________________________

UPDATE SURFACE CONDITION DURING WIND EROSION