# Sediment Modeling

Rivers contribute 95% of sediment entering the ocean. A river’s sediment load is related to many watershed-based processes generally consisting of erosion and deposition. Previous global sediment studies use coarse resolution (0.5 degree) spatial data and empirical models based on area, relief and temperature to derive estimates. Using global datasets of soils, rainfall, elevation and vegetation at a 30 arc-second (approximately 1km) spatial resolution, we were able to use the RUSLE model to predict a global erosion potential map. By summing the erosion potential by basin, we derived estimates of annual sediment load for each river. The estimates for 121 rivers were compared to observed sediment loads and were found to predict within 2 orders of magnitude in all but 3 cases. Techniques were explored to adjust the estimates for the effects of dams, reservoirs and other depositional factors but these did not increase the predictive power of the model.

## Methods

### Overview

The Universal Soil Loss
Equation was developed by the USDA Agricultural Research Service in
1965 as an empirical method to compute soil loss in agricultural
fields. The Revised Universal Soil Loss Equation
(RUSLE) is based on the same empirical theory as USLE but has been
improved in many areas such as methods for calculating rainfall and
slope parameters(**Brooks, 2003**). The RUSLE model is detachment capacity limited
(**Mitasova, **); it only takes into account erosional factors and does
not consider the effects of sediment deposition. Annual erosion is
calculate by combining spatially explicit raster maps according to this
formula:

**E = R K LS C P**

Where:

E = Annual Erosion

R = Rainfall Erosivity Factor

K = Soil Erodibility Factor

LS = Slope Factor

C= Vegetative Cover

P = Management Factor

### Preprocessing Data

All data referred to in subsequent sections have been preprocessed in order to normalize the spatial referencing. Because the model relies on area calculations, it was essential to use an equal-area projection. Seven regions (roughly corresponding to the continents) were defined and each uses a custom Lambert Azimuthal Equal Area projection to minimize distortion. Thus global model is run in seven steps and requires each input dataset to be reprojected into each custom projection and clipped to the region's extents. In order to minimize resampling errors, each region uses a different cell size; the approximate size of a 30 arc second cell at the center lattitude of the projection.

### Slope Factor

The SRTM30 Dataset (****) was used to calculate the slope factor (LS). Since the DEM was clipped to the sea mask and slope calculations require a 3x3 window, it was necessary to add a 5 pixel buffer with a value of zero around the existing DEM to prevent edge effects. Flow accumulation was calculated using the D-infinity algorithm. Slope percentages were calculated using a 3x3 moving window derivative approach (****). LS was calculated from slope and flow accumulation according to the modified method for complex terrain (Mitasova et. al. 1996, 1999)

### Vegetative Management Factor

The MODIS land cover dataset (**http://duckwater.bu.edu/lc/mod12q1.html**), classified according to the IGBP scheme, was used to calculate the vegetative management factor (CP). Each class was assigned a CP factor; certain values were known ( ) and others were interpolated by their relative influence on erosion. The IGBP raster was then reclassified according to these values to derive a CP map.

IGBP Vegetation Class | Vegetative Management Factor (VM) |
---|---|

Water | 0.000 |

Permanent Wetlands | 0.000 |

Snow and Ice | 0.000 |

Evergreen Needleleaf Forest | 0.005 |

Evergreen Broadleaf Forest | 0.007 |

Mixed Forests | 0.009 |

Deciduous Needleleaf Forest | 0.010 |

Deciduous Broadleaf Forest | 0.012 |

Closed Shrublands | 0.050 |

Woody Savannas | 0.050 |

Open Shrublands | 0.075 |

Grasslands | 0.075 |

Savannas | 0.080 |

Urban and Built-Up | 0.100 |

Cropland/Natural Vegetation Mosaic | 0.140 |

Barren or Sparsely Vegetated | 0.200 |

Cropland | 0.320 |

### Soil Erodibility

K was calculated based on
the USDA global soil suborder dataset, available through the NRCS (**CITE**). Each soil order was assigned a factor
based on their described soil properties (**CITE
http://soils.usda.gov/technical/classification/taxonomy/** ) using the
soil nomograph (**USDA, 1980**). The soil suborders were collapsed into orders and reclassified by their K factor.

Soil Order | K |
---|---|

Ocean | 0.00 |

Shifting Sand | 0.99 |

Rock | 0.01 |

Ice | 0.00 |

Gelisols | 0.30 |

Histosols | 0.20 |

Spodosols | 0.20 |

Andisols | 0.30 |

Oxisols | 0.13 |

Verisols | 0.15 |

Aridisols | 0.32 |

Ultisols | 0.18 |

Mollisols | 0.38 |

Alfisols | 0.13 |

Inceptisols | 0.36 |

Entisols | 0.16 |

### Rainfall Erosivity

Using interpolated monthly rainfall data (***) spanning 40 years, the average rainfall erosivity index (R) was calculated. The modified fourniers index (MFI) has shown to be a good predictor of R. The annual MFI for each of the 40 years was calculated by

sum( montly precipitation / total annual precipiatation) over each month

Raster algebra was used on the annual MFI grids to derive an average annual MFI. Finally the grid was interpolated to the optimal regional resolution using regularized spline tension (Mitasova, 1993).

### Validation

Using the Database of World Rivers and their Sediment Yields (**cite**), we calculated the total sediment load (Yield x Area) for each river system. These rivers were observed over varying time spans and the observed values represent long term average sediment loading. In cases where multiple observation stations existed on a single river, the station furthest downstream was used to represent the sediment load at the coast. We joined the observed sediment data with the gtn basins (**cite**) yielding a vector polygon dataset of 123 rivers with known sediment values.

R used to correlate, calculate deviation (modvobs)

### Attempts to improve sediment model by incorporating auxillary data

Using the esri global dams dataset, we performed a spatial intersection with our validation basins and the total dam volume per basin was calculated. Neither the total dam volume nor the volume-to-area ratio were significantly correlated with our models deviations. The Global Lakes and Waterbodies database (GLWD) was also used to calculate the total volume of lakes and resevoirs per basin. This dataset also showed no significant correlation with our model deviation. Despite the fact that dams and resevoirs have a known effect on sediment trapping, neither global dataset had a significant effect on the model.

An additional problem is that sediment transport in rivers is a highly variable phenomenon with possible deviations of 2 orders of magnitude between years. Therefore, an annual average sediment load approach is not likely to predict the sediment loading for a given year.

## Citations

Rueda, C.A., Greenberg, J.A., and Ustin, S.L. StarSpan: A Tool for Fast Selective Pixel Extraction from Remotely Sensed Data. (2005). Center for Spatial Technologies and Remote Sensing (CSTARS), University of California at Davis, Davis, CA.

Mitasova, H., Mitas, L., 1999. Modeling soil detachment with RUSLE 3d using GIS. http://skagit.meas.ncsu.edu/~helena/gmslab/erosion/usle.html. Helena Mitasova, Lubos Mitas, University of Illinois at Urbana-Champaign.

Mitasova, H., J. Hofierka, M. Zlocha, L.R. Iverson, 1996, Modeling topographic potential for erosion and deposition using GIS. Int. Journal of Geographical Information Science, 10(5), 629-641. (reply to a comment to this paper appears in 1997 in Int. Journal of Geographical Information Science, Vol. 11, No. 6)Mitasova, H. and Mitas, L. (1993) Interpolation by Regularized Spline with Tension: I. Theory and Implementation, Mathematical Geology, 25, 641-655.