Erosion Module
Description:
Simulates daily soil erosion, and (optionally) the effect of soil loss on the soil
profile. Code and this document derived from PERFECT.
Methods
The calculation of daily soil loss is performed by one of two submodels:
- A modified MUSLE (Freebairn and Wockner) cover concentration function determined from
QDPI field data to predict soil movement from the inter-contour bank area for clay soils in situations
where peak discharge cannot be adequately predicted. It accounts for variation in soil loss with cover
and runoff volume (the main factors that can be managed), and uses the MUSLE slope-length, erodibility
and practise factors to provide generality. The model has the following form:
Where
E Event soil loss (t/ha)
COV Cover (%)
Q Event runoff (mm)
K Soil erodibility factor (t/ha/EI 30 )
LS Slope length and steepness factor
P Supporting practice factor
The LS factor from Wishmeier and Smith (1978) is related to catchment
slope and length by:
Where
LS Slope length and steepness factor
S Sine of slope angle
L Length of catchment (m)
- A simplified version of the sediment concentration function contained in GUESS
(Carroll). The equation is given as:
Where
E Event soil loss (t/ha)
S Sine of slope angle
cover Fractional surface cover (0-1)
Q Event runoff (mm)
Factor approximating efficiency of entrainment
The cover term cover in the rose equation does not fully account for the
effect of cover. Therefore, the efficiency of entrainment
is further modified by cover by Rose (1985):
Where
Factor approximating efficiency of entrainment
Efficiency of entrainment (bare surface)
COV Surface cover (%)
The submodels parameters are
(rose_lambda), and the ‘-0.15' exponential term (rose_b2).
Feedback from erosion on the soil profile
The module is capable of eroding the soil profile as soil loss occurs (see
profile_reduction). The module remains ignorant of what other modules may be doing through the profile,
all it does is send a delta to whichever module owns dlayer - typically soilwat2. It is the
responsibility of other modules to ensure they remain consistent with this new profile.
The calculation of "dlt_depth_mm" describing the change in each profile layer
for the given amount of soil loss is given by:
dlt_depth_mm(i) = (100.0*soil_loss) / (1000.0*bd(i))
Where
| soil_loss |
Event soil loss (t/ha) |
| bd |
Bulk density of the respective profile layer (g/cc) |
The bulk density for each layer is used to maintain mass balance of soil movement up
through the (constant) soil layer thicknesses.
Erosion continues replacing the lowest profile layer with “imaginary” material of the
same properties until a limit (see bed_depth) is reached. Once the lowest profile layer rests against
bedrock, further soil erosion reduces the thickness of this layer until it is a fraction (see
profile_layer_merge) of its original thickness. Once this fraction is reached, the layer is merged
into the next highest layer. (Currently, many apsim modules cannot sensibly deal with the last
scenario.)
This process continues until the parameter minimum_depth is reached, at which the
simulation stops with a fatal error.
Splitting soil loss into suspended and bed loads.
There is provision to split soil loss into two components: suspended and bed loads. This
reflects how soil loss is measured. Each soil loss submodel is performed twice: once for suspended
ediment, and once for heavy particles.
For the modified MUSLE model, the parameter K is replaced by K_bed and K_susp
( NB. It's important to remove or comment out the original K parameter, otherwise the
module assumes you are still trying to use a single soil loss equation)
For the rose model, parameters lambda and b2 are replaced by lambda_bed, lambda_susp,
b2_bed and b2_susp. Once again, you must remove the original definitions of lambda and b2.
Reporting the separate losses is through soil_loss_bed, soil_loss_susp, both in units of
t/ha, and sed_conc_bed, sed_conc_susp (g/l). soil_loss returns total soil loss, and sed_conc returns
total sediment concentration.
Communicating with water balance modules
As the erosion module was derived substantially from PERFECT (ie. a daily timestep), it
does not use the advanced features of APSIM's apswim and surface modules (eg. hydrographs or peak runoff
rate). SWIM users should be aware that peak runoff rate is not used in calculation of soil erosion at
this time. Surface users should be aware that surface roughness factors do not enter erosion or
runoff.
The calculation of cover used in the soil loss equation should be the same cover as used
in runoff prediction. When soilwat provides water balance, this cover is exported to the system (and
thus erosion) as “total_cover”.
However, SWIM does not provide this variable. If the surface module is active, it will
provide a variable “surf_cover”, which is the combined crop and residue covers on the soil surface. To
use this cover as the cover term in the soil loss equation, add the following to the manager's rules:
xxx.manager.init
total_cover = 0.0
xxx.manager.start_of_day
total_cover = surface.surf_cover
To use erosion with apswim, add the following rules:
xxx.manager.init
total_cover = 0.0
xxx.manager.start_of_day
total_cover = 1.0 – (1.0 - crop1.cover_tot) * (1.0 – crop2.cover_tot) * … * (1.0 - residue2.residue_cover)
Procedures:
The module follows these steps at initialisation:
- All variables are set to zero.
- The parameter file is read, if necessary, ls factor is calculated for the freebairn model.
- The static parameters are written to the summary file.
During daily processing, the following procedures occur:
- The daily state variables are reset.
- Today's inputs are requested from the system.
- Soil loss is calculated as a function of cover and runoff.
- If there is soil loss, the profile is changed and sent back to its owner.
- The new depth to bedrock is calculated.
- Control returns to the system.
System interface names
Parameter file:
Name |
Description |
model |
Either “rose” or “freebairn” |
profile_reduction |
Either “on” or “off” |
profile_layer_merge |
Fraction of original size below which the lowest layer is merged into the layer above (0-1) |
minimum_depth |
If the profile erodes to this depth, the simulation is stopped (mm) |
slope |
Slope of plot in percent (%) |
slope_length |
Length of plot (m) |
bed_depth |
depth to bedrock (mm) |
cover_extra (optional) |
fudge factor added to total_cover (-1.0 -1.0) |
Freebairn specific parameters
|
p_factor |
Supporting practise factor (unitless) |
Either |
k_factor |
Soil erodibility factor (t/ha/EI 30 ) |
or |
k_factor_bed |
Soil erodibility factor for bed load (t/ha/EI 30 ) |
k_factor_susp |
Soil erodibility factor for suspended load (t/ha/EI 30 ) |
|
Rose specific parameters
Either |
entrain_eff |
Efficiency of entrainment - bare surface () |
or |
entrain_eff _bed |
Efficiency of entrainment - bare surface - bed load () |
entrain_eff _susp |
Efficiency of entrainment - bare surface - suspended load () |
|
Either |
eros_rose_b2 |
?? |
or |
eros_rose_b2_bed |
?? |
eros_rose_b2_susp |
?? |
|
Inputs from other modules on a daily timestep:
day |
|
year |
|
runoff |
daily runoff (mm) |
total_cover |
combined crop and residue cover (0 - 1) |
bd(mxlayr) |
moist bulk density of soil (g/cm^3) |
dlayer(mxlayr) |
thickness of soil layer i (mm) |
Outputs to other modules as requested:
soil_loss |
todays soil loss (t/ha) |
soil_loss_bed |
todays soil loss in bedload (t/ha) |
soil_loss_susp |
todays soil loss in suspension (t/ha) |
soil_loss_mm |
todays soil loss from the topmost profile layer (mm) |
sed_conc |
todays sediment concentration (g/l) |
sed_conc_bed |
todays sediment concentration in bedload (g/l) |
sed_conc_susp |
todays sediment concentration in suspended load (g/l) |
bed_depth |
todays depth to bedrock (mm) |
erosion_cover |
todays cover used to drive soil loss equations (0 - 1) |
Resets to other modules:
dlt_dlayr() |
thickness of soil layer i (mm) |
Typical, sample parameters
Dataset |
slope |
Slope length |
P factor |
soil erodibility K |
Lambda_bare |
b2 |
|
% |
M |
|
(t/ha/EI30) |
|
|
Greenmount |
6.5 |
60 |
1.0 |
0.38 (1) |
0.77 (2) |
15 (3) |
Greenwood |
4.5 |
37 |
1.0 |
0.38 (1) |
0.70 (3) |
15 (3) |
gmtf.erosion.parameters
model = freebairn (1)
slope = 6.5 (%)
slope_length = 60.0 (m)
bed_depth = 1500. (mm)
profile_reduction = off
profile_layer_merge = 0.05 ()
minimum_depth = 100.0 (mm)
k_factor = 0.38 ()
p_factor = 1.0 ()
gmtr.erosion.parameters
model = rose (2)
slope = 6.5 (%)
slope_length = 60.0 (m)
bed_depth = 1500. (mm)
profile_reduction = off
profile_layer_merge = 0.05 ()
minimum_depth = 100.0(mm)
entrain_eff = 0.77 ()
eros_rose_b2 = 0.15 ()
gwd.erosion.parameters
model = rose (3)
slope = 4.5 (%)
slope_length = 37.0 (m)
bed_depth = 1500. (mm)
profile_reduction = off
profile_layer_merge = 0.05 ()
minimum_depth = 100.0(mm)
entrain_eff = 0.7 ()
eros_rose_b2 = 0.15 ()
(1) Freebairn, Silburn & Loch (1989). Aust.J.Soil Res. 27: 199-211
(2) Silburn & Loch (1992) 5th Aust. Soil con. Conf.
(3) Rose (1985) Adv. in Soil Science Vol 2
Known bugs and defects
Cover isn't resolved. Currently, erosion gets its
cover from “total_cover”, which soilwat calculates as the total crop
and residue cover. This should be runoff_cover, which soilwat doesnt
export at the moment.
Bibliography
Carroll, C., Rose, C.W. and Crawford, S.J. (1986).
GUESS - Theory Manual. School of Australian Environmental Studies,
Griffith University , Queensland ,
Freebairn, D.F. and Wockner, G.H. (1986). A study of soil erosion on Vertisols of the Eastern Darling Downs, Australian Journal of Soil Research , 19 , 133-46
Littleboy, M., Silburn, D.M., Freebairn, D.M.,
Woodruff, D.R., and Hammer, G.L. (1989) PERFECT - A computer simulation
model of Productivity Erosion Runoff Functions to Evaluate Conservation
Techniques. QDPI, Brisbane.
Rose, C.W. (1985). Developments in soil erosion and deposition models. Advances in Soil Science , Volume 2 , 1-63.
Wischmeier, W.H. and Smith, D.D. (1978). Predicting rainfall erosion losses, a guide to conservation planning. Agricultural Handbook number 537 , USDA, 58pp.
Programmers notes
The calculation of soil loss ends up in two
variables, soil_loss_bed and soil_loss_susp. When the user describes a
single soil loss equation, the result is left in soil_loss_bed, with
soil_loss_susp left at 0. All reporting is done as the sum of these two
variables.
Related Material
Nutrient reduction:
Nutrient reduction is represented by "moving" the
soil layers downwards through the soil profile (the same as moving the
nutrient pools of the soil layers upwards through the soil). Layers
take on the nutrient characteristics of the layer below in proportion
to the depth increment gained from that layer. If there is no
"bedrock", all layers would eventually all end up with the nutrient
properties of the subsoil (bottom layer). If the accumulated depth loss
due to erosion exceeds the depth between the bottom of the profile and
the depth to bedrock (and profile reduction is on), the bottom layer
does not contain nutrients from below and its nutrient content will
diminish (as nutrients are exported upwards through the soil)
The algorithm for "top down" removal is:
1. determine the loss for each layer. For the top layer, this is
given by
conc_kg_kg = variable(1)/(1000.0*bd(1)*dlayr(1)*10.0)
loss_kg_ha = soil_loss * 1000.0* enr * conc_kg_kg
where
conc_kg_kg = kg of nutrient per kg of soil (kg/kg)
variable = content of particular nutrient (kg/ha),
bd = bulk density in (g/cm^3)
dlayr = depth of layer in (mm)
soil_loss = soil loss (t/ha)
enr = enrichment ratio.
"enr" is calculated from
enr = enr_a_coeff * (1000.0 * soil_loss)**(-1.0 * enr_b_coeff)
and is bounded to (1.0 < enr < enr_a_coeff):
enr = amin1(enr_a_coeff, enr)
enr = amax1(enr, 1.0)
For the remaining layers, the loss from each layer is that which moved into the layer above;
layer_loss = variable(i) *
layer_divide(dlt_depth_mm(i), dlayr(i))
where i is the index of the layer we are working on,
and layer_divide(a,b) bounds the result of (a/b) to 0.0 -> 1.0. This
prevents taking away more than is present in the layer.
"dlt_depth_mm(i)" is an array of deltas describing the change in each
profile depth for the given amount of soil loss, given by:
dlt_depth_mm(i) = (100.0*soil_loss) / (1000.0*bd(i))
The bulk density for each layer is used to maintain
mass balance of soil movement up through the (constant) soil layer
thicknesses.
The gain to each layer is found in a similar manner, ie, the fraction of the layer below that moves into this layer:
layer_gain = variable(i+1) *
layer_divide(dlt_depth_mm(i+1), dlayr(i+1))
Obviously, the lowest layer cannot use this formula,
as there is no layer beneath, and we have to cater for bedrock. The
procedure is to assume that the lowest layer gains a fraction of
itself, so long as it is not resting on bedrock:
layer_gain = variable(i) *
layer_divide(dlt_depth_mm(i), dlayr(i))
! check we're not going into bedrock
if (suml(dlayr, mxlayr) +
dlt_depth_mm(n_layrs) .gt. bed_depth .and.
profile_reduction .eq. on ) then
layer_gain = 0.0
At this point, we find the new value of the variable for this layer in the new profile:
variable(i) = variable(i) + layer_gain - layer_loss
Two more steps follow,:-
1. If profile reduction is on, check to see if the bottom layer is too thin, and merge it into the next layer up if so.
2. Perform a mass balance check:
\Sum{yesterdays variable} + loss from layer 1 - gain to lowest layer
= \Sum{Todays variable}
This procedure is repeated for all the nitrogen / carbon variables:
snh4(mxlayr)
inert_c(mxlayr)
biom_c(mxlayr)
biom_n(mxlayr)
hum_c(mxlayr)
hum_n(mxlayr)
fom_n(mxlayr)
fpool1(mxlayr)
fpool2(mxlayr)
fpool3(mxlayr)
The total amounts of N and C lost on eroded soil is calculated by summing losses from the relevant N and C pools:
n_loss_in_sed = snh4_loss + biom_n_loss +
hum_n_loss + fom_n_loss
c_loss_in_sed = biom_c_loss + hum_c_loss +
fpool1_loss + fpool2_loss + fpool3_loss
----
Finally, reducing the profile layer depth (dlayr)
takes place following a similar method to before; by reducing from the
bottom until the profile rests against bedrock:
tot_depth = suml(dlayr, n_layrs) + dlt_depth_mm(n_layrs)
if (tot_depth .gt. bed_depth ) then
overrun = tot_depth - bed_depth
do 2000 i = n_layrs, 1, -1
if (overrun .gt. 0.0) then
if(overrun .le. dlayr(i)) then
dlayr(i) = dlayr(i) - overrun
overrun = 0.0
else
overrun = overrun - dlayr(i)
dlayr(i) = 0.0
endif
endif
2000 continue
endif
and the lowest profile layer is merged if necessary.
The threshold for merging layers is specified in the
parameter file as a proportion of the original layer, so after a merge,
this threshold must be recalculated as a proportion of the next upper
layer.