Erosion

Simulates daily soil erosion, and (optionally) the effect of soil loss on the soil profile. Code and this document are derived from PERFECT.

Methods:

The calculation of daily soil loss is performed by one of two submodels:

Sub Model 1: Freebairn

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:

{Equation 1}

 

Where;

• E: Event soil loss (t/ha)
• COV: Cover (%)
• Q: Event runoff (mm)
• K: Soil erodibility factor (t/ha/EI30 )
• 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:

{Equation 2}

Where;

• LS: Slope length and steepness factor
• S: Sine of slope angle
• L: Length of catchment (m)

Dimensional Correctness of Modified MUSLE

Note that the equation for E is not dimensionally correct (the units of E (t/ha) do not match those on the right hand side of the equation).  To be dimensionally correct E should be in units of (mm t)/(ha EI30) instead of t/ha. The existance of the EI30 unit can be traced back to the original USLE(Wischmeier and Smith, 1978) equation on which this equation is based.  In the USLE EI30 unit is cancelled out because both R (Soil Erosivity) and K (Soil Erodibility) have EI30 as units. The unit remains in the new equation because E (Erosion) is calculated using Q (Runoff) not R (Soil Erosivity). The new equation still uses K however, which means the EI30 unit from the K is no longer cancelled out. For consistency, K is kept in the same units in the new equation. so that users who have measured the values of K for use in USLE can use the same measured values (in the same units) in the new equation. However, keeping these same units for consistency results in the new equation being dimensionally incorrect.  The other extra dimension, the mm unit, is also gained as a consequence of using Q (Runoff) instead of R.

To resolve these dimensional issues we firstly define K as a being unit less however with the same value as for K used in USLE.

From first principles the equation above starts out as,

 

In USLE equation it is R x K which gives us Erosion in t/ha.

Therefore we need to replace R such that (R_runoff x Q) x K still equals t/ha of erosion.

K is now unitless however,

Therefore R_runoff must be in t/(ha x mm) which means R_runoff x Q(mm) x K = t/ha

R_runoff is therefore the tons of soil lost per ha per millimeter of runoff.

 

E = R_runoff x Q x K x LS x P

OR

E = (sediment erosion in t/(ha x mm)) x Q x K x LS x P

 

From R_runoff, we can remove a unit conversion for (grams of water as runoff) to (mm of water as runoff)

and reduce R_runoff to a fraction of (grams of soil lost(sediment)) per (gram of water as runoff).

 

sediment concentration = 1 gram of sediment / 1 gram of water

 

1 gram of sediment = (1/10^6) tons of sediment

 

1000Kg of Water = 1m^3

10^6 grams of water = 1000mm x 1m^2

1gram of water = (1000 x (mm x m^2))/10^6

 

1 ha = 100m x 100m = 10^4m^2

1m^2 = 1/(10^4)

 

1 gram of water = (1000 x (mm x ha)) / (10^6 x 10^4)

1 gram of water = (1/10^7) x (mm x ha)

 

frac_to_(t/ha mm) =  1 gram of sediment/ 1 gram of water = (tonnes/(10^6)) x (10^7 / (mm x ha) = 10 t/(mm x ha)

 

We are now left with

 

E = sediment concentration function x (10 t/(ha mm)) x Q(mm) x K x LS x P

 

 

 

The remaining sediment concentration function however is going to be a function of the cover fraction.

Rather than using cover fraction in the function it would be nicer to use the cover percentage and

just move this conversion out of the function so that the function now has cover as a percentage.

 

E =  [sediment concentration function x percent_to_frac x frac_to_(t/(ha mm))] x Q(mm) x K x LS x P

 

E  = [sediment concentration function x (1/100) x (10 t/(ha mm))] x Q(mm) x K x LS x P

 

E  = [sediment concentration function x (1/10) t/ha/mm] x Q(mm) x K x LS x P

 

This matches the simpler form at the start of this section {Equation 1}
where,

sediment concentration function = 16.52 – 0.46 * cover + 0.0031 cover^2

and (1/10) actually has its units removed,

(1/10) t/(ha mm) = (1/10)

We now see Equation 1 in its full, dimensionally correct form:

E (t/ha) = 16.52 – 0.46 * cover + 0.0031 cover^2 x (1/10)(t/(ha mm)) x LS x K x P x Q(mm)

 

(when cover is less than 50%)

 

AND

 

E (t/ha) = -0.0254 * cover + 2.54 x (1/10)(t/(ha mm)) x LS x K x P x Q(mm)

 

(when cover is greater than 50%)

QDPI field data was used to obtain the form and the coefficients for sediment concentration function. This is possible because E, Q, K, LS, P are all known via measurement.

As mentioned above, We do not publish the full dimensionally correct form.
For consistency sake we replace the unitless K with the same dimensions as the USLE definition of K.
We then ignore the t/(ha mm) units on the (1/10) as well.
This results in a dimensionally incorrect form of the equation but one that is easier to understand for users of the USLE.

Right Units for K

In the section above about dimensional correctness we stated,
“To resolve these dimensional issues we firstly define K as a being unitless however with the same value as for K used in USLE.”

Even though when using the modified MUSLE (Freebairn and Wockner) we don’t care about the units of K
(because to keep modified MUSLE (Freebairn and Wockner) dimensionally correct we drop off the units)
the value of K still needs to be the value in the correct units.

As a result of the history of the USLE there are 3 possible units for K.
Originally the USLE equation was in Imperial units. Or rather the units of R and K were in Imperial units. (because all the other variables in USLE are unitless)
The Equation was later converted to Metric units. Or rather the units for R and K were converted to Metric units.
Finally the Equation was converted to SI units (which are metric but are in the Official SI metric units). So once again R and K were converted into SI units.

At this point we need to state that the correct units of K for the modified MUSLE (Freebairn and Wockner) is K in Metric units.

So if you have K in either Imperial units or SI units,YOU MUST CONVERT IT TO METRIC BEFORE USING IT.
Note: There is no R in the modified MUSLE (Freebairn and Wockner) so we don’t need to worry about converting the units of R, like we would if we were using the USLE.

To convert K(SI) to K(Metric) simply multiply the K(SI) by Gravity (9.81 m/s^2).

The sensible range of values for K(SI)
is between 0.001 and 0.09

The sensible range of values for K(Metric)
is between 0.01 and 0.88

 

 

 

You can derive the above conversion from SI to Metric as follows,

According to the USLE,
the units of K (the Erodability of the Soil) are in t/(ha EI30)
where EI30 is the unit of R (the Erosivity of the Soil).
The idea of this being that the units of R which are EI30 will be cancelled out by the EI30 in K and just leave the t/ha for the result of the equation (because all the other variables in USLE are unitless).

Therefore if we wish to convert K(SI) into K(Metric)
we first must work out the conversion of units of EI30(SI) to EI30(Metric) (or to state another way, we must first work out the conversion of units of R(SI) to R(Metric))
because EI30(SI) are part of the units of K(SI).

Since,

R(SI) = EI30(SI) = (MJ/ha) x (mm/h)

R(Metric) = EI30(Metric) = 100m x (t/ha) x (cm/h)

 

EI30 is a measure of Kinetic Energy of falling rain and Rainfall Intensity.

Joules is the unit of the Kinetic Energy of the falling rain.

We divide this Energy over the area it falls which is why we get MJ/ha.

The mm/h is the Maximum Rainfall intensity falling over 30min period.

 

Joules = Work = Force x Distance =  (Mass x Acceleration) x Distance = Mass x Gravity x Distance = Kg x 9.81 (m/s^2) x m

 

Starting with R(SI) we want to get to R(Metric)

R(SI) = EI30(SI) = (MJ/ha) x (mm/h)

 

Therefore,

 

using

Joules = Kg x 9.81 (m/s^2) x m

 

 

R(SI) = 10^6 x Kg x 9.81 (m/s^2) x m x (1/ha) x (mm/h)

 

R(SI)/(9.81 m/s^2) = 10^6 x Kg x m x (1/ha) x (mm/h)

 

= 10^3 x t x m x (1/ha) x (mm/h)

= 10^2 x m x (t/ha) x (10 x mm/h)

= 100m x (t/ha) x (cm/h)

= R(Metric)

 

hence,

R(SI)/ (9.81 m/s^2)  = R(Metric)

or

R(SI)/Gravity = R(Metric)

 

 

or

EI30(SI)/Gravity = EI30(Metric)

 

 

Now that we have worked out the conversion of R(SI) to R(Metric) we can now work out the conversion of K(SI) to K(Metric).

K(SI) = t/(ha x EI30(SI))

K(Metric) = t/(ha x EI30(Metric))

 

using, EI30(Metric) = EI30(SI)/Gravity

 

K(Metric) = t/(ha (EI30(SI)/Gravity) )

 

K(Metric) = (t x Gravity)/(ha x EI30(SI))

 

K(Metric) = Gravity x t/(ha x EI30(SI))

 

K(Metric) = Gravity x K(SI)

 

K(Metric) = 9.81 (m/s^2) x K(SI)

 

Sub Model 2: Rose

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:

1. All variables are set to zero.
2. The parameter file is read, if necessary, ls factor is calculated for the freebairn model.
3. The static parameters are written to the summary file.

During daily processing, the following procedures occur:

1. The daily state variables are reset.
2. Today’s inputs are requested from the system.
3. Soil loss is calculated as a function of cover and runoff.
4. If there is soil loss, the profile is changed and sent back to its owner.
5. The new depth to bedrock is calculated.
6. Control returns to the system.

System interface names

Name	Description	Units
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	(%)
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 (model = ‘freebairn’)

Name	Description	Units
p_factor	Supporting practise factor	(unitless)

 

Either

Name	Description	Units
k_factor	Soil erodibility factor	(t/ha/EI 30 )

or

Name	Description	Units
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 (model = ‘rose’)

Either

Name	Description	Units
entrain_eff	Efficiency of entrainment – bare surface	()

or

Name		Description	Units
entrain_eff _bed		Efficiency of entrainment – bare surface – bed load	()
entrain_eff _susp	Efficiency of entrainment – bare surface – suspended load	()

Either

Name	Description
eros_rose_b2	??

or

Name	Description
eros_rose_b2_bed	??
eros_rose_b2_susp	??

 

Inputs from other modules on a daily timestep:

Name	Description	Units
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:

Name	Description	Units
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:

Name	Description	Units
dlt_dlayr()	thickness of soil layer i	(mm)

dlt_dlayr is delta layer depth or change in layer depth.

It allows you to alter the layer depths from their current values by a certain amounts specified in this array variable.

ie. if

dlayer = 10 20 30 40 50

and

dlt_dlayr = 6 -7 -8 9 10

then dlayer would become,

dlayer = 16 13 22 49 60

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)*

* See the corresponding examples below.

The following are examples of Apsim Erosion parameters in the Published Scientific Literature:
(1) Freebairn, Silburn & Loch (1989). Aust.J.Soil Res. 27: 199-211

[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 ()

(2) Silburn & Loch (1992) 5th Aust. Soil con. Conf.

[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 ()

(3) Rose (1985) Adv. in Soil Science Vol 2

[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 ()

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.
• Freebairn, D.M., Loch, R.J., Silburn, D.M., 1996. Chapter 9 Soil erosion and soil conservation for vertisols, in: Ahmad, N., Mermut, A. (Eds.), Developments in Soil Science, Vertisols and Technologies for Their Management. Elsevier, pp. 303–362.

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.

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 

1. Determine the loss for each layer.

1.1 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)

1.2 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.

2. The gain to each layer is found in a similar manner, ie, the fraction of the layer below that moves into this layer:

2.2

layer_gain = variable(i+1) *  layer_divide(dlt_depth_mm(i+1), dlayr(i+1))

2.3 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

3. At this point, we find the new value of the variable for this layer in the new profile:

3.1

variable(i) = variable(i) + layer_gain – layer_loss

Two more steps follow,:-

3.2 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.

3.3 Perform a mass balance check:

\Sum{yesterdays variable} + loss from layer 1 – gain to lowest layer  = \Sum{Todays variable}

3.4 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

 

4. 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.