Garden with Insight v1.0 Help: Erosion - Wind Erosion

The original EPIC wind erosion model (WEQ) required daily mean wind speed as a driving variable. The new EPIC wind erosion model WECS (Wind Erosion Continuous Simulation) requires the daily distribution of wind speed to take advantage of the more mechanistic erosion equation. The new approach estimates potential wind erosion for a smooth bare soil by integrating the erosion equation through a day using the wind speed distribution. Potential erosion is adjusted using four factors based on soil properties, surface roughness, cover, and distance across the field in the wind direction.

Basic Equation

The basic WECS wind erosion equation is [Equation 139] where YW is the wind erosion in kg/m, FI is the soil erodibility factor, FR is the surface roughness factor, FV is the vegetative cover factor, FD is the mean unsheltered travel distance of wind across a field, DW is the duration of wind greater than the threshold velocity in sec, and YWR is the wind erosion rate in kg/m*sec at time t.

Equation 139

YW = FI * FR * FV * FD * (integration from 0 to DW of) YWR
Code:
YW = 8640.0 * FI * FR * FV * ROKF * (integration from 0 to DW of) (YWR * FD)
note that FD is multiplied into the integrated YWR (see equation 141)
also code bounds this aResult at maxWindErosionPerDay_tPha, an input
Variables:
YW = WindErosion_tPha
FI = soilErodibilityFactorWE
FR = surfaceRoughnessFactorWE
FV = vegetativeCoverFactorWE
FD = meanUnsheltDistanceWE
DW = durationOfWindOverThresholdVelocityWE_sec
YWR = windErosionRateForFractionOfDay_kgPmsec
ROKF = coarseFragmentFactorWE

The wind erosion rate is calculated with the equation of Skidmore(1986)

YWR = C * (rho(a) / g) * (v*(o)2 - v*(tau)2 - 0.5 * (SW/WP)2)3/2 (Equation 140)

where C is a parameter ~~2.5, rho(a) is the air density in kg/m3, g is the acceleration of gravity in m/s2, v*(o) is the friction velocity in m/sec, and v*(tau) is the threshold friction velocity in m/sec. SW and WP are the actual and 1500 kPa water contents of the top soil layer (10 mm thick), respectively. The soil water term of the equation was developed by Chepil (1956) and Skidmore (1986). Substituting acceleration of gravity (9.8 m/sec2) and assuming air density is 1 kg/m3 gives [Equation 141].

Equation 141

YWR = 0.255 * power(sqr(v*(o)) - sqr(v*(tau)) - 0.5 * sqr(SW/WP), 3/2)
Code:
ustw = sqr(v*(tau)) + 0.5 * sqr(SW/WP) + 0.5 - tlmf
YWR = 0.255 * power(sqr(v*(o)) - ustw, 3/2)
YWR = 0.255 * power(sqr(v*(o)) - (sqr(v*(tau)) + 0.5 * sqr(SW/WP) + 0.5 - tlmf), 3/2)
YWR = 0.255 * power(sqr(v*(o)) - sqr(v*(tau)) - 0.5 * sqr(SW/WP) - 0.5 + tlmf, 3/2)
so the code adds a -0.5 and a tlmf term to the part inside the power
code also adds meanUnsheltDistanceFactorWE multiplier (FD), but does not include it
in the final equation (so it works out, see equation 139)
code also adds check: aResult is 0 if v*(o) - ustw < 0
Variables:
YWR = WindErosionRateForFractionOfDay_kgPmsec
v*(o) = frictionVelocityWE_mPsec
ustw = weForFractionOfDayFactor
FD = meanUnsheltDistanceFactorWE

The friction velocity is estimated with the equation

v*(o) = K * v / (log(Z/Z(o)) (Equation 142)

where K is Von Karmans constant ~~0.4, v is the wind speed in m/sec at height Z (10 m for WECS input), and Z(o) is the aerodynamic roughness (0.00055 m used in WECS). Substituting the WECS values into equation 142 gives [Equation 143].

Equation 143

v(o) = 0.0408 * v
Code:
same
Variables:
v*(o) = FrictionVelocityWE_mPsec
v = windSpeedForFractionOfDay_mPsec

The threshold friction velocity is estimated with the equation

v*(tau) = 0.1 * sqrt(g * (delta-p/rho(a)) * D) (Equation 144)

where delta-p is the mineral soil density (2650 kg/m3) and D is the soil particle size in m. Substituting constants into equation 144 and expressing D in micrometers gives [Equation 145].

Equation 145

v*(tau) = 0.0161 * sqrt(D)
Code:
same
Variables:
v*(tau) = ThresholdFrictionVelocityWE_mPsec
D = soilParticleDiameter_microns

Soil Erodibility Factor

WECS uses the soil erodibility concept of WEQ expressed in dimensionless form with the equation [Equation 146] where I is the soil erodibility factor of the Woodruff and Siddoway (1965) model in t/ha and FI is the dimensionless soil erodibility factor of the new model.

Equation 146

FI = I / 695.0
Code:
based on soil sand, silt and clay
above equation may be the same. do not know how I is derived.
Variables:
FI = SoilErodibilityFactorWE
I = soilErodibilityFactorWEWoodruffAndSiddoway_tPha

Roughness Factor

The surface roughness factor (FR) is based upon the shelter angle developed by Potter et al. (1990). This roughness index calculates the erodible fraction of the soil surface by estimating the portion susceptible to abrasion by saltating particles. The shelter angle index incorporates both roughness due to random cloddiness and oriented roughness (ridges) due to tillage operations. The effect of oriented roughness varies as a function of wind direction, which is selected each day so that the statistical distribution of wind direction approaches that of the simulation site. FR is estimated with the equation [Equation 147] where wn(1) is the descent angle of saltating sand grains (about 15 degrees from horizontal). A 15 degree impact angle has also been shown to cause maximum aggregate abrasion (Hagen et al., 1988).

Equation 147

FR = 1.0 - exp(-power(wn(1) / RFB, RFC))
wn(1) = 15 degrees
Code:
same except code uses 5.0 for wn(1)
Variables:
FR = SurfaceRoughnessFactorWE
wn(1) = 5.0
RFB = coeffSurfaceRoughnessFactorWE
RFC = expSurfaceRoughnessFactorWE

The coefficient RFC is calculated with the equation [Equation 148] where RHTT is the ridge height in mm.

Equation 148

RFC = 0.77 * power(1.002, RHTT)
Code:
same
Variables:
RFC = ExpSurfaceRoughnessFactorWE
RHTT = ridgeHeight_mm

The coefficient RFB is estimated with the equations [Equation 149], [Equation 150] and [Equation 151] where RRF is the clod roughness factor, RRUF is the random roughness in mm, RIF is the ridge roughness factor, and wn(2) is the angle of the wind relative to ridges. Both RRUF and RHTT are altered by wind and water erosion and tillage.

Equation 149

RFB = RRF + RIF
Code:
same except for lower bound at 1.0
Variables:
RFB = CoeffSurfaceRoughnessFactorWE
RRF = clodRoughnessFactorWE
RIF = ridgeRoughnessFactorWE

Equation 150

RRF = 11.9 * (1.0 - exp(-power(RRUF / 9.8, 1.3)))
Code:
same
Variables:
RRF = ClodRoughnessFactorWE
RRUF = randomRoughnessWE_mm

Equation 151

RIF = abs(sin(wn(2) * 1.27 * power(RHTT, 0.52))
Code:
Variables:
RIF = RidgeRoughnessFactorWE
wn(2) = angleOfWindRelativeToRidges_rad
RHTT = ridgeHeight_mm

Vegetative Cover Factor

The vegetative cover factor is based on the approach used in the original EPIC model. A vegetative cover equivalent factor is simulated daily as a class function of standing live biomass, standing dead residue, and flat crop residue [Equation 152] where VE is the vegetative cover equivalent factor, SB is the standing biomass in t/ha, SR is the standing crop residue in t/ha, FR is the flat residue in t/ha, and w(1), w(2) and w(3) are crop specific coefficients.

Equation 152

VE = 0.253 * power((w(1) * SB + w(2) * SR + w(3) * FR), 1.36)
Code:
VE = w(1) * SB + w(2) * SR + w(3) * FR
Variables:
VE = VegCoverEquivFactorWE
SB = standingLiveBiomass_tPha
SR = standingDeadResidue_tPha
FR = surfaceLayerFlatCropResidue_tPha
w(1) = windErosionFactorStandingLive
w(2) = windErosionFactorStandingDead
w(3) = windErosionFactorFlatResidue

The vegetative equivalent is converted to a vegetative factor for use in the new model [Equation 153]. Equation 153 assures that 0.0 <= FV <= 1.0.

Equation 153

FV = VE / (VE + exp(0.48 - 1.32 * VE))
Code:
FV = 1.0 - VE / (VE + exp(0.48 - 1.32 * VE))
Variables:
FV = VegetativeCoverFactorWE_frn
VE = vegCoverEquivFactorWE

Unsheltered Distance Factor

Field length along the prevailing wind direction is calculated as in the original model (Cole at al., 1982) by considering the field dimensions and orientation and the wind direction. [Equation 154] where WL is the unsheltered field length along the prevailing wind direction in m, FL is the field length in m, FW is the field width in m, theta is the wind direction clockwise from north in radians, and phi is the clockwise angle between field length and north in radians.

Equation 154

WL = FL * FW / (FL * cos(pi/2 + theta - phi) + FW * sin(pi/2 + theta - phi))
Code:
same
Variables:
WL = UnsheltFieldLengthAlongPrevailingWindDir_m
FL = fieldLength_m
FW = fieldWidth_m
theta = windDirectionForDay_rad
phi = fieldLengthOrientationFromNorth_rad

The new model distance factor (FD) is calculated as described by Stout (1990) using the equation [Equation 155] where wn(3) is a parameter determined experimentally to lie in the range 50.0 < wn(3) < 90.0. A value of 70.0 is used in EPIC.

Equation 155

FD = 1.0 - exp(-WL / wn(3))
wn(3) = 70.0
Code:
same except numerator is 0.07, not 70
Variables:
FD = MeanUnsheltDistanceFactorWE
WL = unsheltFieldLengthAlongPrevailingWindDir_m
wn(3) = 70.0