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
|