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

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:

Freebairn specific parameters

Rose specific parameters

Inputs from other modules on a daily timestep:

Outputs to other modules as requested:

Resets to other modules:

Typical, sample parameters

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