SDR: Sediment Delivery Ratio¶
Summary¶
The objective of the InVEST Sediment Delivery Ratio (SDR) model is to map overland sediment generation and delivery to the stream. In the context of global change, such information can be used to study the service of sediment retention in a catchment. This is of particular interest for reservoir management and instream water quality, both of which may be economically valued.
Introduction¶
Erosion and overland sediment retention are natural processes that govern the sediment concentration in streams. Sediment dynamics at the catchment scale are mainly determined by climate (in particular rain intensity), soil properties, topography, and vegetation; and anthropogenic factors such as agricultural activities or dam construction and operation. Main sediment sources include overland erosion (soil particles detached and transported by rain and overland flow), gullies (channels that concentrate flow), bank erosion, and mass erosion (or landslides; see Merrit 2003 for a review). Sinks include on-slope, floodplain or instream deposition, and reservoir retention, as summarized in Figure 1. Conversion of land use and changes in land management practices may dramatically modify the amount of sediment running off a catchment. The magnitude of this effect is primarily governed by: i) the main sediment sources (land use change will have a smaller effect in catchments where sediments are not primarily coming from overland flow); and ii) the spatial distribution of sediment sources and sinks (for example, land use change will have a smaller effect if the sediment sources are buffered by vegetation).
Increases in sediment yield are observed in many places in the world, dramatically affecting water quality and reservoir management (UNESCO 2009). The sediment retention service provided by natural landscapes is of great interest to water managers. Understanding where the sediments are produced and delivered allow managers to design improved strategies for reducing sediment loads. Changes in sediment load can have impacts on downstream irrigation, water treatment, recreation and reservoir performance.
Outputs from the sediment model include the sediment load delivered to the stream at an annual time scale, as well as the amount of sediment eroded in the catchment and retained by vegetation and topographic features. Note that SDR only creates biophysical results. For valuation of the sediment retention service, appropriate valuation approaches will be highly dependent on the particular application and context, and need to be implemented independently of InVEST.
Fig. 9 General catchment sediment budget. The relative size of the arrows changes depending on the environment. The InVEST model focuses on the overland sources and sinks, and does not include the others.¶
The Model¶
Sediment Delivery¶
The sediment delivery module is a spatially-explicit model working at the spatial resolution of the input digital elevation model (DEM) raster. For each pixel, the model first computes the amount of annual soil loss from that pixel, then computes the sediment delivery ratio (SDR), which is the proportion of soil loss actually reaching the stream. Once sediment reaches the stream, we assume that it ends up at the catchment outlet, thus no in-stream processes are modeled. This approach was proposed by Borselli et al. (2008) and has received increasing interest in recent years (Cavalli et al., 2013; López-vicente et al., 2013; Sougnez et al., 2011). See the User Guide section Differences between the InVEST SDR model and the original approach developed by Borselli et al. (2008) for further discussion.
Annual Soil Loss¶
The amount of annual soil loss on pixel \(i\), \(usle_i\) (units: \(tons\cdot ha^{-1} yr^{-1}\)), is given by the revised universal soil loss equation (RUSLE1):
where
\(R_i\) is rainfall erosivity (units: \(MJ\cdot mm (ha\cdot hr\cdot yr)^{-1})\),
\(K_i\) is soil erodibility (units: \(ton\cdot ha\cdot hr (MJ\cdot ha\cdot mm)^{-1}\)),
\(LS_i\) is a slope length-gradient factor (unitless)
\(C_i\) is a cover-management factor (unitless)
and \(P_i\) is a support practice factor (Renard et al., 1997). (cf. also in (Bhattarai and Dutta, 2006)). (unitless)
The \(LS_i\) factor is given from the method developed by Desmet and Govers (1996) for a two-dimensional surface:
where
\(S_i\) is the slope factor for grid cell \(i\) calculated as a function of slope radians \(\theta\)
\[\begin{split}S = \left\{\begin{array}{lr} 10.8\cdot\sin(\theta)+0.03, & \text{where } \theta < 9\% \\ 16.8\cdot\sin(\theta)-0.50, & \text{where } \theta \geq 9\% \\ \end{array}\right\}\end{split}\]\(A_{i-in}\) is the contributing area (\(m^2\)) at the inlet of a grid cell which is computed from the Multiple-Flow Direction method
\(D\) is the grid cell linear dimension (\(m\))
\(x_i\) is the mean of aspect weighted by proportional outflow from grid cell \(i\) determined by a Multiple-Flow Direction algorithm. It is calculated by
\[x_i = \sum_{d\in{\{0,7\}}} x_d\cdot P_i(d)\]where \(x_d = |\sin \alpha(d)| + |\cos \alpha(d)|\), \(\alpha(d)\) is the radian angle for direction \(d\), and \(P_i(d)\) is the proportion of total outflow at cell \(i\) in direction \(d\).
\(m\) is the RUSLE length exponent factor.
To avoid overestimation of the LS factor in heterogeneous landscapes, long slope lengths are capped to a maximum value of 122m that is adjustable as a user parameter (Desmet and Govers, 1996; Renard et al., 1997).
The value of \(m\), the length exponent of the LS factor, is based on the classical USLE, as discussed in (Oliveira et al., 2013):
Sediment Delivery Ratio¶
Step 1. Based on the work by Borselli et al. (2008), the model first computes the connectivity index (\(IC\)) for each pixel. The connectivity index describes the hydrological linkage between sources of sediment (from the landscape) and sinks (like streams.) Higher values of \(IC\) indicate that source erosion is more likely to make it to a sink (i.e. is more connected), which happens, for example, when there is sparse vegetation or higher slope. Lower values of \(IC\) (i.e. lower connectivity) are associated with more vegetated areas and lower slopes.
\(IC\) is a function of both the area upslope of each pixel (\(D_{up}\)) and the flow path between the pixel and the nearest stream (\(D_{dn}\)). If the upslope area is large, has lower slope, and good vegetative cover (so a low USLE C factor), \(D_{up}\) will be low, indicating a lower potential for sediment to make it to the stream. Similarly, if the downslope path between the pixel and the stream is long, has lower slope and good vegetative cover, \(D_{dn}\) will be low.
\(IC\) is calculated as follows:
Figure 2. Conceptual approach used in the model. The sediment delivery ratio (SDR) for each pixel is a function of the upslope area and downslope flow path (Equations 3, 4, 5).
Thresholded slopes \(S_{th}\) and cover-management factors \(C_{th}\) are used in calculating \(D_{up}\) and \(D_{dn}\). A lower bound is set to avoid infinite values for \(IC\). An upper bound is also applied to the slope to limit bias due to very high values of \(IC\) on steep slopes. (Cavalli et al., 2013).
\(D_{up}\) is the upslope component defined as:
where \(\bar{C}_{th}\) is the average thresholded \(C\) factor of the upslope contributing area, \(\bar{S}_{th}\) is the average thresholded slope gradient of the upslope contributing area (\(m/m\)) and \(A\) is the upslope contributing area (\(m^2\)). The upslope contributing area is delineated from a Multiple-Flow Direction algorithm.
The downslope component \(D_{dn}\) is given by:
where \(d_i\) is the length of the flow path along the ith cell according to the steepest downslope direction (\(m\)) (see Figure 2), \(C_{th, i}\) and \(S_{th, i}\) are the thresholded cover-management factor and the thresholded slope gradient of the ith cell, respectively. Again, the downslope flow path is determined from a Multiple-Flow Direction algorithm.
Step 2. The SDR ratio for a pixel \(i\) is then derived from the conductivity index \(IC\) following (Vigiak et al., 2012):
where \(SDR_{max}\) is the maximum theoretical SDR, set to an average value of 0.8 (Vigiak et al., 2012), and \(IC_0\) and \(k\) are calibration parameters that define the shape of the SDR-IC relationship (which is an increasing function). The effect of \(IC_0\) and \(k\) on the SDR is illustrated below:
Figure 3. Relationship between the connectivity index IC and the SDR. The maximum value of SDR is set to \(SDR_{max}=0.8\). The effect of the calibration are illustrated by setting \(k_b=1\) and \(k_b=2\) (solid and dashed line, respectively), and \(IC_0=0.5\) and \(IC_0=2\) (black and grey dashed lines, respectively).
Bare Soil¶
The SDR for bare soil is is calculated in the same way, simply leaving out the cover-management factor:
where \(IC_{bare}\) is the connectivity index for bare soil, defined as:
\(D_{up, bare}\) is the upslope component for bare soil, defined as:
and \(D_{dn, bare}\) is the downslope component for bare soil, defined as:
Sediment Export¶
The sediment export from a given pixel i \(E_i\) (units: \(tons\cdot ha^{-1} yr^{-1}\)), is the amount of sediment eroded from that pixel that actually reaches the stream. Sediment export is given by:
The total catchment sediment export \(E\) (units: \(ton\cdot ha^{-1} yr^{-1}\)) is given by:
\(E\) is the value used for calibration/validation purposes, in combination with other sediment sources, if data are available.
Sediment Downslope Deposition¶
This model also makes an estimate of the amount of sediment that is deposited on the landscape downstream from the source that does not reach the stream. Knowing the spatial distribution of this quantity will allow users to track net change of sediment on a pixel (gain or loss) which can inform land degradation indices.
Sediment export to stream from pixel \(i\) is defined in equation (56). The other component of the mass balance from the USLE is that sediment which does not reach the stream. This sediment load must be deposited somewhere on the landscape along the flowpath to the stream and is defined as follows
Due to the nature of the calculation of SDR, the quantity \(E_i\) has accounted for the downstream flow path and biophysical properties that filter sediment to stream. Thus, we can model the flow of \(E'_i\) downstream independently of the flow of \(E_i\).
To do this, we assume the following properties about how \(E_i\) and SDR behave across a landscape:
Property A: SDR monotonically increases along a downhill flowpath: As a flowpath is traced downhill, the value of SDR will monotonically increase since amount of downstream flow distance decreases. Note there is the numerical possibility that a downstream pixel has the same SDR value as an upstream pixel. The implication in this case is that no on-pixel sediment flux deposition occurs along that step.
Property B: All non-exporting sediment flux on a boundary stream pixel is retained by that pixel: If pixel \(i\) drains directly to the stream there is no opportunity for further downstream filtering of \(E_i\). Since \(E_i\) is the inverse of Ei, the implication is that the upstream flux (defined as Fi below) must have been deposited on the pixel.
Given these two properties, we see that the amount of \(E_i\) retained on a pixel must be a function of:
the absolute difference in SDR values from pixel \(i\) to the downstream pixel(s) drain, and
how numerically close the downstream SDR value is to 1.0 (the stream pixel).
These mechanics can be captured as a linear interpolation of the difference of pixel i’s SDR value with its downstream SDR counterpart with respect to the difference of pixel i’s difference with a theoretical maximum downstream SDR value 1.0. Formally,
The \(d\) in \(dR_i\) indicates a delta difference and \(p(i,k)\) is the proportion of flow from pixel \(i\) to pixel \(j\). This notation is meant to invoke the intution of a derivative of \(Ri\). Note the boundary conditions are satisfied:
In the case of Property A (downstream \(SDR_k=SDR_i\)), the value of \(dR_i=0\) indicating no \(F_i\) will be retained on the pixel.
In the case of Property B (downstream \(SDR_k=1\) because it is a stream) the value of \(dR_i=1\) indicating the remaining \(F_i\) is retained on the pixel.
Now we define the amount of sediment flux that is retained on any pixel in the flowpath using \(dR_i\) as a weighted flow of upstream flux:
where \(F_j\) is the amount of sediment-export-that-does-not-reach-the stream “flux”, defined as:
Streams and Optional Drainage Layer¶
The model’s stream map is the union of the calculated stream layer and the input drainage layer (if provided). The model calculates a stream layer (stream.tif) by thresholding the flow accumulation raster (flow_accumulation.tif) by the threshold flow accumulation (TFA) value:
(62)¶\[\begin{split}stream_{TFA,i} = \left\{\begin{array}{lr} 1, & \text{if } flow\_accum_{i} \geq TFA \\ 0, & \text{otherwise} \\ \end{array}\right\}\end{split}\]
If the optional drainage input is provided, the model includes it (stream_and_drainage.tif):
(63)¶\[stream_{drainage,i} = stream_{TFA,i} \text{ OR } stream_{input,i}\]
The final stream layer (\(stream_{TFA}\), or \(stream_{drainage}\) if the optional drainage input is provided) is used to determine \(d_i\) for the SDR calculations.
In some situations, the index of connectivity defined by topography does not represent actual flow paths, which may be influenced by artificial connectivity instead. For example, sediments in urban areas or near roads are likely to be conveyed to the stream with little retention. The (optional) drainage raster identifies the pixels that are artificially connected to the stream, irrespective of their geographic position (e.g. their distance to the stream network). Pixels from the drainage layer are treated similarly to pixels of the stream network; in other words, the downstream flow path will stop at pixels of the drainage layer (and the corresponding sediment load will be added to the total sediment export).
Defined Area of Outputs¶
SDR and several other model outputs are defined in terms of distance to stream (\(d_i\)). Therefore, these outputs are only defined for pixels that drain to a stream on the map (and so are within the streams’ watershed). Pixels that do not drain to any stream will have nodata in these outputs. The affected output files are: d_dn.tif, ic.tif, e_prime.tif, sdr_factor.tif, sdr_bare_soil.tif, d_dn_bare_soil.tif, ic_bare_soil.tif, sed_retention.tif. sed_retention_index.tif, sediment_deposition.tif, and sed_export.tif
If you see areas of nodata in these outputs that can’t be explained by missing data in the inputs, it is likely because they are not hydrologically connected to a stream on the map. This may happen if your DEM has pits or errors, if the map boundaries do not extend far enough to include streams in that watershed, or if your threshold flow accumulation value is too high to recognize the streams. Check the stream output (stream.tif) and make sure that it aligns as closely as possible with the streams in the real world.
Example: Below is an example of the effect of threshold flow accumulation on the defined extent, in an area with multiple watersheds that are not hydrologically connected. The top row shows streams (stream.tif), while the bottom row shows SDR (sdr_factor.tif).
In the left column, with a TFA value of 100, streams exist in both the bottom-left and top-right watersheds. The SDR raster is defined everywhere that the inputs are defined except for a small patch on the right edge that does not drain to any stream.
In the right column, with a TFA value of 1000, there are no streams at all in the upper-right watershed. As a result, pixels in that watershed do not drain to any stream, and the corresponding SDR raster is undefined in that area.
Limitations¶
Among the main limitations of the model is its reliance on the USLE (Renard et al., 1997). This equation is widely used but is limited in scope, only representing rill/inter-rill erosion processes. Other sources of sediment include gully erosion, streambank erosion, and mass erosion. A good description of the gully and streambank erosion processes is provided by Wilkinson et al. 2014, with possible modeling approaches. Mass erosion (landslide) is not represented in the model but can be a significant source in some areas or under certain land use change, such as road construction.
A corollary is that the descriptions of the impact on ecosystem services (and any subsequent valuation) should account for the relative proportion of the sediment source from the model compared to the total sediment budget (see the section on Evaluating sediment retention services).
In addition, as an empirical equation developed in the United States, the USLE has shown limited performance in other areas – even when focusing on sheet and rill erosion. Based on local knowledge, users may modify the soil loss equation implemented in the model by altering the R, K, C, P inputs to reflect findings from local studies (Sougnez et al., 2011).
The model is very sensitive to the k and IC0 parameters, which are not physically based. The emerging literature on the modeling approach used in the InVEST model (Cavalli et al., 2013; López-vicente et al., 2013; Sougnez et al., 2011; Vigiak et al., 2012) provides guidance to set these parameters, but users should be aware of this limitation when interpreting the model’s absolute values.
Given the simplicity of the model and low number of parameters, outputs are very sensitive to most input parameters. Errors in the empirical parameters of the USLE equations will therefore have a large effect on predictions. Sensitivity analyses are recommended to investigate how the confidence intervals in input parameters affect the study conclusions.
Differences between the InVEST SDR model and the original approach developed by Borselli et al. (2008)¶
The InVEST SDR model is based on the concept of hydrological connectivity, as parameterized by Borselli et al. (2012). This approach was selected since it requires a minimal number of parameters, uses globally available data, and is spatially explicit. In a comparative study, Vigiak et al. (2012) suggested that the approach provides: “(i) large improvement in predicting specific sediment yields, (ii) ease of implementation, (iii) scale-independency; and (iv) a formulation capable of accounting for landscape variables and topology in line with sedimentological connectivity concepts”. The approach has also been used to predict the effect of land use change (Jamshidi et al., 2013).
The following points summarize the differences between InVEST and the Borselli model:
The weighting factor is directly implemented as the USLE C factor (other researchers have used a different formulation, e.g. roughness index based on a high-resolution DEM (Cavalli et al., 2013))
The \(SDR_{max}\) parameter used by Borselli et al. is set to 0.8 by default to reduce the number of parameters. Vigiak et al. (2012) propose to define \(SDR_{max}\) as the fraction of topsoil particles finer than coarse sand (<1 mm).
Evaluating Sediment Retention Services¶
Sediment Retention Services¶
Translating the biophysical impacts of altered sediment delivery to human well-being metrics depends very much on the decision context. Soil erosion, suspended sediment and deposited sediment can have both negative and positive impacts on various users in a watershed (Keeler et al, 2012). These include, but are not limited to:
Reduced soil fertility and reduced water and nutrient holding capacity, impacting farmers
Increase in treatment costs for drinking water supply
Reduced lake clarity diminishing the value of recreation
Increase in total suspended solids impacting health and distribution of aquatic populations
Increase in reservoir sedimentation diminishing reservoir performance or increasing sediment control costs
Increase in harbor sedimentation requiring dredging to preserve harbor function
Sediment Retention¶
Sediment retention is computed as follows:
Sediment Retention Index¶
An index of sediment retention is computed by the model as follows:
which represents the avoided soil loss by the current land use compared to bare soil, weighted by the SDR factor. This index underestimates retention since it does not account for the retention from upstream sediment flowing through the given pixel. Therefore, this index should not be interpreted quantitatively. We also note that in some situations, index values may be counter-intuitive: for example, urban pixels may have a higher index than forest pixels if they are highly connected to the stream. In other terms, the SDR (second factor) can be high for these pixels, compensating for a lower service of avoided soil loss (the first factor): this suggests that the urban environment is already providing a service of reduced soil loss compared to an area of bare soil.
Quantitative Valuation¶
An important note about assigning a monetary value to any service is that valuation should only be done on model outputs that have been calibrated and validated. Otherwise, it is unknown how well the model is representing the area of interest, which may lead to misrepresentation of the exact value. If the model has not been calibrated, only relative results should be used (such as an increase of 10%) not absolute values (such as 1,523 tons, or 42,900 dollars.)
Sediment retention at the subwatershed level From a valuation standpoint, an important metric is the difference in retention or yield across scenarios. For quantitative assessment of the retention service, the model uses as a benchmark a hypothetical scenario where all land is cleared to bare soil: the value of the retention service is then based on the difference between the sediment export from this bare soil catchment and that of the scenario of interest. This output is termed “sed_retent” in the watershed summary table and sed_retention.tif in the raster outputs. Similarly, the sediment retention provided by different user-provided scenarios may be compared with the baseline condition (or each other) by taking the difference in sediment export between scenario and baseline. This change in export can represent the change in sediment retention service due to the possible future reflected in the scenario.
Additional sources and sinks of sediment As noted in the model limitations, the omission of some sources and sinks of sediment (gully erosion, stream bank erosion, and mass erosion) should be considered in the valuation analyses. In some systems, these other sources of sediment may dominate and large changes in overland erosion may not make a difference to overall sediment concentrations in streams. In other words, if the sediment yields from two scenarios differ by 50%, and the part of rill/inter-rill erosion in the sediment budget in 60%, then the actual change in erosion that should be valued for avoided reservoir sedimentation is 30% (50% x .6).
One complication when calculating the total sediment budget is that changes in climate or land use result in changes in peak flow during rain events, and are thus likely to affect the magnitude of gully and streambank erosion. While the magnitude of the change in other sediment sources is highly contextual, it is likely to be in the same direction as the change in overland erosion: a higher sediment overland transport is indeed often associated with higher flows, which likely increase gully and bank erosion. Therefore, when comparing across scenarios, the absolute change may serve as a lower bound on the total impact of a particular climate or land use change.
Appendix 2 summarizes options to represent the additional sources and sinks of erosion in the model.
Replacement and avoided cost frameworks, versus willingness to pay approaches With many ecosystem service impacts, and sediment impacts in particular, monetary valuation is relatively simple if an avoided mitigation cost or replacement cost method is deemed appropriate. In this situation, beneficiaries are assumed to incur a cost that is a function of the biophysical metric (e.g., suspended sediment increases treatment costs). However, it is important to recognize that the avoided cost or replacement cost approaches assume the mitigating actions are worthwhile for the actor undertaking them. For example, if a reservoir operator deems that the costs associated with dredging deposited sediment are not worth the benefits of regaining lost storage capacity, it is not appropriate to value all deposited sediment at the unit cost of dredging. Similarly, an increase in suspended sediment for drinking water supplies may be met by increasing treatment inputs or switching to an alternate treatment technology. Avoiding these extra costs could then be counted as economic benefits. However, in some contexts, private water users may decide that the increase in sediment content is acceptable, rather than incur additional treatment expenses. They are economically worse off, but by not paying for additional treatment, the replacement cost approach becomes an upper bound on their economic loss. Their economic loss is also no longer captured by their change in financial expenditures, which further complicates the analysis.
Note, however, that this bounding approach may be entirely appropriate for initial assessment of the significance of different benefit streams, i.e. if the most expensive approach does not have a significant impact, then there is no need to refine the analysis to utilize more detailed approaches such as willingness-to-pay (for consumers) or impacts on net revenues (for producers). However, if the impact is large and there is no good reason to believe that the relevant actors will undertake the mitigating activities, then a willingness-to-pay framework is the appropriate path to take. For an introduction to the techniques available, see http://ecosystemvaluation.org/dollar_based.htm.
Time considerations Generally, economic and financial analysis will utilize some form of discounting that recognizes the time value of money, benefits, and use of resources. Benefits and costs that accrue in the future “count for less” than benefits and costs that are borne close to the present. It is important that any economic or financial analysis be cognizant of the fact that the SDR model represents only average annual impacts under steady state conditions. This has two implications for valuation. First, users must recognize that the impacts being valued may take some time to come about: It is not the case that the full steady state benefits would begin accruing immediately, even though many of the costs might. Second, the annual averaging means that cost or benefit functions displaying nonlinearities on shorter timescales should (if possible) be transformed, or the InVEST output should be paired with other statistical analysis to represent important intra- or inter-annual variability.
Data Needs¶
Raster inputs may have different cell sizes, and they will be resampled to match the cell size of the DEM. Therefore, all model results will have the same cell size as the DEM.
Workspace (required). Folder where model outputs will be written. Make sure that there is ample disk space, and write permissions are correct.
Suffix (optional). Text string that will be appended to the end of output file names, as “_Suffix”. Use a Suffix to differentiate model runs, for example by providing a short name for each scenario. If a Suffix is not provided, or changed between model runs, the tool will overwrite previous results.
Digital elevation model (DEM) (required). Raster dataset with an elevation value for each cell. Make sure the DEM is corrected by filling in sinks, and compare the output stream maps with hydrographic maps of the area. To ensure proper flow routing, the DEM should extend beyond the watersheds of interest, rather than being clipped to the watershed edge. [units: meters]
Rainfall erosivity index (R) (required). Raster dataset, with an erosivity index value for each cell. This variable depends on the intensity and duration of rainfall in the area of interest. The greater the intensity and duration of the rain storm, the higher the erosion potential. The erosivity index is widely used, but in case of its absence, there are methods and equations to help generate a grid using climatic data. [units: \(MJ\cdot mm\cdot (ha\cdot h\cdot yr)^{-1}\)]
Soil erodibility (K) (required). Raster dataset, with a soil erodibility value for each cell. Soil erodibility, K, is a measure of the susceptibility of soil particles to detachment and transport by rainfall and runoff. [units: \(tons\cdot ha\cdot h\cdot (ha\cdot MJ\cdot mm)^{-1}\)]
Land use/land cover (LULC) (required). Raster dataset, with an integer LULC code for each cell. All values in this raster MUST have corresponding entries in the Biophysical table.
Watersheds (required). A shapefile of polygons. This is a layer of watersheds such that each watershed contributes to a point of interest where water quality will be analyzed. Format: An integer field named ws_id is required, with a unique integer value for each watershed.
Biophysical table (required). A .csv (Comma Separated Value) table containing model information corresponding to each of the land use classes in the LULC raster. All LULC classes in the LULC raster MUST have corresponding values in this table. Each row is a land use/land cover class and columns must be named and defined as follows:
lucode: Unique integer for each LULC class (e.g., 1 for forest, 3 for grassland, etc.) Every value in the LULC map MUST have a corresponding lucode value in the biophysical table.
usle_c: Cover-management factor for the USLE, a floating point value between 0 and 1.
usle_p: Support practice factor for the USLE, a floating point value between 0 and 1.
Threshold flow accumulation (required). The number of upstream cells that must flow into a cell before it is considered part of a stream, which is used to classify streams from the DEM. This threshold directly affects the expression of hydrologic connectivity and the sediment export result: when a flow path reaches the stream, sediment deposition stops and the sediment exported is assumed to reach the catchment outlet. It is important to choose this value carefully, so modeled streams come as close to reality as possible. See Appendix 1 for more information on choosing this value. Integer value, with no commas or periods - for example “1000”.
\(k_b\) and \(IC_0\) (required): Two calibration parameters that determine the shape of the relationship between hydrologic connectivity (the degree of connection from patches of land to the stream) and the sediment delivery ratio (percentage of soil loss that actually reaches the stream; cf. Figure 3). The default values are \(k_b=2\) and \(IC_0=0.5\).
\(\mathbf{SDR_{max}}\) (required): The maximum SDR that a pixel can reach, which is a function of the soil texture. More specifically, it is defined as the fraction of topsoil particles finer than coarse sand (1000 μm; Vigiak et al. 2012). This parameter can be used for calibration in advanced studies. Its default value is 0.8.
\(\mathbf{l_{max}}\) (required): The maximum allowed value of the L parameter when calculating the LS factor. Calculated values that exceed this are clamped to this value. Its default value is 122 but reasonable values in literature place it anywhere between 122-333 see Desmet and Govers, 1996 and Renard et al., 1997.
Drainage layer (optional) A raster with 0s and 1s, where 1s correspond to pixels artificially connected to the stream (by roads, stormwater pipes, etc.) and 0s are assigned to all other pixels. The flow routing will stop at these “artificially connected” pixels, before reaching the stream network, and the corresponding sediment exported is assumed to reach the catchment outlet.
Interpreting Results¶
The resolution of the output rasters will be the same as the resolution of the DEM provided as input.
[Workspace] folder:
Parameter log: Each time the model is run, a text (.txt) file will be created in the Workspace. This file will list the parameter values and output messages for that run and will be named according to the service, the date and time, and the suffix. When contacting NatCap about errors in a model run, please include the parameter log.
rkls.tif (type: raster; units: tons/pixel): Total potential soil loss per pixel in the original land cover from the RKLS equation. Equivalent to the soil loss for bare soil. (Eq. (44), without applying the \(C\) or \(P\) factors)
sed_export.tif (type: raster; units: tons/pixel): The total amount of sediment exported from each pixel that reaches the stream. (Eq. (56))
sediment_deposition.tif (type: raster; units: tons/pixel): The total amount of sediment deposited on the pixel from upstream sources as a result of retention. (Eq. (60))
stream_and_drainage.tif (type: raster): If a drainage layer is provided, this raster is the union of that layer with the calculated stream layer(Eq. (63)). Values of 1 represent streams, values of 0 are non-stream pixels. Compare this layer with a real-world stream map, and adjust the Threshold Flow Accumulation so that this map matches real-world streams as closely as possible.
usle.tif (type: raster; units: tons/pixel): Total potential soil loss per pixel in the original land cover calculated from the USLE equation. (Eq. (44))
sed_retention.tif (type: raster; units: tons/pixel): Map of sediment retention with reference to a watershed where all LULC types are converted to bare ground. (Eq. (64))
sed_retention_index.tif (type: raster; units: tons/pixel, but should be interpreted as relative values, not absolute): Index of sediment retention, used to identified areas contributing more to retention with reference to a watershed where all LULC types are converted to bare ground. This is NOT the sediment retained on each pixel (see Section on the index in “Evaluating Sediment Retention Services” above). (Eq. (65))
watershed_results_sdr.shp: Table containing biophysical values for each watershed, with fields as follows:
sed_export (units: tons/watershed): Total amount of sediment exported to the stream per watershed. This should be compared to any observed sediment loading at the outlet of the watershed. Knowledge of the hydrologic regime in the watershed and the contribution of the sheetwash yield into total sediment yield help adjust and calibrate this model. (Eq. (57) with sum calculated over the watershed area)
usle_tot (units: tons/watershed): Total amount of potential soil loss in each watershed calculated by the USLE equation. (Sum of USLE from (44) over the watershed area)
sed_retent (units: tons/watershed): Difference in the amount of sediment delivered by the current watershed and a hypothetical watershed where all land use types have been converted to bare ground. (Sum of (64) over the watershed area)
sed_dep (units: tons/watershed): Total amount of sediment deposited on the landscape in each watershed, which does not enter the stream. (Sum of \(R_i\) from (60) over the watershed area)
[Workspace]\intermediate_outputs folder:
cp.tif: \(C\cdot P\) factor (Eq. (44)), derived by mapping usle_c and usle_p from the biophysical table to the LULC raster.
d_dn_bare_soil.tif: downslope factor of the index of connectivity, ignoring the cover-management factor as if the soil were bare (Eq. (55))
d_dn.tif: downslope factor of the index of connectivity (Eq. (50))
d_up_bare_soil.tif: upslope factor of the index of connectivity, ignoring the cover-management factor as if the soil were bare (Eq. (54))
d_up.tif: upslope factor of the index of connectivity (Eq. (49))
e_prime.tif: sediment downslope deposition, the amount of sediment from a given pixel that does not reach a stream (Eq. (58))
f.tif: sediment flux for sediment that does not reach the stream (Eq. (61))
stream.tif: stream raster calculated directly from flow accumulation, flow direction, and the TFA value (Eq. (62)).
flow_accumulation.tif: flow accumulation, derived from flow direction
flow_direction.tif: MFD flow direction. Note: the pixel values should not be interpreted directly. Each 32-bit number consists of 8 4-bit numbers. Each 4-bit number represents the proportion of flow into one of the eight neighboring pixels.
ic_bare_soil.tif: index of connectivity, ignoring the cover-management factor as if the soil were bare (Eq. (53))
ic.tif: index of connectivity (Eq. (46))
ls.tif: LS factor for USLE (Eq. (45))
pit_filled_dem.tif: DEM after any pits are filled
s_accumulation.tif: flow accumulation weighted by the thresholded slope. Used in calculating s_bar.
s_bar.tif: mean thresholded slope gradient of the upslope contributing area (\(\bar{S}_{th}\) in eq. (49))
s_inverse.tif: inverse of the thresholded slope (\(1/S_{th}\) in eq. (50))
sdr_bare_soil.tif: sediment delivery ratio, ignoring the cover-management factor as if the soil were bare (Eq. (52))
sdr_factor.tif: sediment delivery ratio (Eq. (51))
slope.tif: slope in radians, calculated from the pit-filled DEM
slope_threshold.tif: slope in radians, thresholded to be no less than 0.005 and no greather than 1 (eq. (47))
w_threshold.tif: cover-management factor thresholded to be no less than 0.001 (eq. (48))
w_accumulation.tif: flow accumulation weighted by the thresholded cover-management factor. Used in calculating w_bar.
w_bar.tif: mean thresholded cover-management factor for upslope contributing area (\(\bar{C}_{th}\) in eq. (49))
w.tif: cover-management factor derived by mapping usle_c from the biophysical table to the LULC raster
weighted_avg_aspect.tif: average aspect weighted by flow direction (\(x\) in eq. (45))
ws_inverse.tif: Inverse of the thresholded cover-management factor times the thresholded slope (\(1/(C_{th} \cdot S_{th})\) in eq. (50))
Comparison with Observations¶
The sediment yield (sed_export) predicted by the model can be compared with available observations. These can take the form of sediment accumulation in a reservoir or time series of Total Suspended Solids (TSS) or turbidity. In the former case, the units are the same as in the InVEST model (tons per year). For time series, concentration data need to be converted to annual loads (LOADEST and FLUX32 are two software facilitating this conversion). Time series of sediment loading used for model validation should span over a reasonably long period (preferably at least 10 years) to attenuate the effect of inter-annual variability. Time series should also be relatively complete throughout a year (without significant seasonal data gaps) to ensure comparison with total annual loads.
A global database of sediment yields for large rivers can be found on the FAO website: http://www.fao.org/nr/water/aquastat/sediment/index.stm Alternatively, for large catchments, global sediment models can be used to estimate the sediment yield. A review of such models was performed by de Vente et al. (2013).
A key thing to remember when comparing modeled results to observations is that the model represents rill-inter-rill erosion only. As indicated in the Introduction three other sources of sediment may contribute to the sediment budget: gully erosion, stream bank erosion, and mass erosion. The relative importance of these processes in a given landscape needs to be determined to ensure appropriate model interpretation.
For more detailed information on comparing with observations, and associated calibration, see Hamel et al (2015).
If there are dams on streams in the analysis area, it is possible that they are retaining sediment, such that it will not arrive at the outlet of the study area. In this case, it may be useful to adjust for this retention when comparing model results with observed data. For an example of how this was done for a study in the northeast U.S., see Griffin et al 2020. The dam retention methodology is described in the paper’s Appendix, and requires knowing the sediment trapping efficiency of the dam(s).
Appendix 1: Data Sources¶
Digital Elevation Model¶
Land Use/Land Cover¶
Watersheds¶
Threshold Flow Accumulation¶
Rainfall Erosivity Index (R)¶
R should be obtained from published values, as calculation is very tedious. For calculation, R equals the annual average of EI values, where E is the kinetic energy of rainfall (in \(MJ\cdot ha^{-1}\)) and I30 is the maximum intensity of rain in 30 minutes (in mm.hr-1). A review of relationships between precipitation and erosivity index around the world is provided by Renard and Freimund (1994).
General guidance to calculate the R index can be found in the FAO Soils bulletin 70 (Roose, 1996): http://www.fao.org/3/t1765e/t1765e0e.htm. It is also possible that area- or country-specific equations for R have been derived, so it is worth doing a literature search for these.
A global map of rainfall erosivity (30 arc-seconds, ~1km at the equator) is available from the European Commission: https://esdac.jrc.ec.europa.eu/content/global-rainfall-erosivity.
In the United States, national maps of the erosivity index can be found through the United States Department of Agriculture (USDA) and Environmental Protection Agency (EPA) websites. The USDA published a soil loss handbook (https://www3.epa.gov/npdes/pubs/ruslech2.pdf ) that contains a hard copy map of the erosivity index for each region. Using these maps requires creating a new line feature class in GIS and converting to raster. Please note that conversion of units is also required: multiplication by 17.02 is needed to convert from US customary units to MJ.mm.(ha.h.yr)-1, as detailed in Appendix A of the USDA RUSLE handbook (Renard et al., 1997).
The EPA has created a digital map that is available at https://archive.epa.gov/esd/archive-nerl-esd1/web/html/wemap_mm_sl_rusle_r_qt.html. The map is in a shapefile format that needs to be converted to raster, along with an adjustment in units.
Soil Erodibility (K)¶
Texture is the principal factor affecting K, but soil profile, organic matter and permeability also contribute. It varies from 70/100 for the most fragile soil to 1/100 for the most stable soil (in US customary units). Erodibility is typically measured on bare reference plots, 22.2 m-long on 9% slopes, tilled in the direction of the slope and having received no organic matter for three years.
Global soil data are available from the Soil and Terrain Database (SOTER) Programme (https://data.isric.org:443/geonetwork/srv/eng/catalog.search). They provide some area-specific soil databases, as well as SoilGrids globally.
The FAO also provides global soil data in their Harmonized World Soil Database: https://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/, but it is rather coarse.
In the United States free soil data is available from the U.S. Department of Agriculture’s NRCS gSSURGO, SSURGO and gNATSGO databases: https://www.nrcs.usda.gov/wps/portal/nrcs/main/soils/survey/geo/. They also provide ArcGIS tools (Soil Data Viewer for SSURGO and Soil Data Development Toolbox for gNATSGO) that help with processing these databases into spatial data that can be used by the model. The Soil Data Development Toolbox is easiest to use, and highly recommended if you use ArcGIS and need to process U.S. soil data.
Please note that conversion of units may be required: multiplication by 0.1317 is needed to convert from US customary units to \(ton\cdot ha\cdot hr\cdot (ha\cdot MJ\cdot mm)^{-1}\), as detailed in Appendix A of the USDA RUSLE handbook (Renard et al., 1997).
Alternatively, the following equation can be used to calculate K (Renard et al., 1997):
In which K = soil erodibility factor (\(t\cdot ha\cdot hr\cdot (MJ\cdot mm\cdot ha)^{-1}\); M = (silt (%) + very fine sand (%))(100-clay (%)) a = organic matter (%) b = structure code: (1) very structured or particulate, (2) fairly structured, (3) slightly structured and (4) solid c = profile permeability code: (1) rapid, (2) moderate to rapid, (3) moderate, (4) moderate to slow, (5) slow and (6) very slow.
When profile permeability and structure are not available, soil erodibility can be estimated based on soil texture and organic matter content, based on the work of Wischmeier, Johnson and Cross (reported in Roose, 1996). The OMAFRA fact sheet summarize these values in the following table (http://www.omafra.gov.on.ca/english/engineer/facts/12-051.htm):
Average OMC |
OMC<2% |
OMC>2% |
|
|---|---|---|---|
Clay |
0.22 |
0.24 |
0.21 |
Clay loam |
0.3 |
0.33 |
0.28 |
Coarse sandy loam |
0.07 |
0.07 |
0.07 |
Fine sand |
0.08 |
0.09 |
0.06 |
Fine sandy loam |
0.18 |
0.22 |
0.17 |
Heavy clay |
0.17 |
0.19 |
0.15 |
Loam |
0.3 |
0.34 |
0.26 |
Loamy fine sand |
0.11 |
0.15 |
0.09 |
Loamy sand |
0.04 |
0.05 |
0.04 |
Loamy very fine sand |
0.39 |
0.44 |
0.25 |
Sand |
0.02 |
0.03 |
0.01 |
Sandy clay loam |
0.2 |
0.2 |
0.2 |
Sandy loam |
0.13 |
0.14 |
0.12 |
Silt loam |
0.38 |
0.41 |
0.37 |
Silty clay |
0.26 |
0.27 |
0.26 |
Silty clay loam |
0.32 |
0.35 |
0.3 |
Very fine sand |
0.43 |
0.46 |
0.37 |
Very fine sandy loam |
0.35 |
0.41 |
0.33 |
The soil erodibility values (K) in this table are in US customary units, and require the 0.1317 conversion mentioned above. Values are based on the OMAFRA Fact sheet. Soil textural classes can be derived from the FAO guidelines for soil description (FAO, 2006, Figure 4).
A special case is the K value for water bodies, for which soil maps may not indicate any soil type. A value of 0 can be used, assuming that no soil loss occurs in water bodies.
Sometimes, soil maps may also have holes in places that are not water bodies (such as glaciers.) Here, look at a land cover map to see what is happening on the landscape. If it is a place where erosion is unlikely to happen (such as rock outcrops), a value of 0 may be used. However, if the area seems like it should have soil data, you can use a nearest neighbor GIS function, or manually set those areas to the dominant soil type that surrounds the missing data.
P and C Coefficients¶
The cover-management factor, C, accounts for the specified crop and management relative to tilled continuous fallow. The support practice factor, P, accounts for the effects of contour plowing, strip-cropping or terracing relative to straight-row farming up and down the slope. These values will need to be obtained from a literature search. Several references on estimating these factors can be found online:
USDA: RUSLE handbook (Renard et al., 1997)
OMAFRA: USLE Fact Sheet http://www.omafra.gov.on.ca/english/engineer/facts/12-051.htm
U.N. Food and Agriculture Organization http://www.fao.org/3/T1765E/t1765e0c.htm
Calibration Parameters \(IC_0\) and \(k_b\)¶
\(IC_0\) and \(k_b\) are calibration parameters that define the relationship between the index of connectivity and the sediment delivery ratio (SDR). Vigiak et al. (2012) suggest that \(IC_0\) is landscape independent and that the model is more sensitive to \(k_b\) . Advances in sediment modeling science should refine our understanding of the hydrologic connectivity and help improve this guidance. In the meantime, following other authors (Jamshidi et al., 2013), we recommend setting these parameters to their default values ( \(IC_0\) =0.5 and \(k_b\) =2), and using \(k_b\) only for calibration (Vigiak et al., 2012).
For more detailed information on sensitivity analysis and calibration, see Hamel et al (2015).
Appendix 2: Representation of Additional Sources and Sinks of Sediment¶
The InVEST model predicts the sediment delivery only from sheetflow erosion, thus neglecting other sources and sinks of sediment (e.g. gully erosion, streambank, landslides, stream deposition, etc.), which can affect the valuation approach. Adding these elements to the sediment budget requires good knowledge of the sediment dynamics of the area and is typically beyond the scope of ecosystem services assessments. General formulations for instream deposition or gully formation are still an area of active research, with modelers systematically recognizing large uncertainties in process representation (Hughes and Prosser, 2003; Wilkinson et al., 2014). Consultation of the local literature to estimate the relative importance of additional sources and sinks is a more practical approach to assess their effect on the valuation approach.
Type of erosion |
Process |
Representation |
|---|---|---|
Gully erosion (connected to the stream) |
If the gully is connected to the stream, it is essentially a small stream tributary and its contribution may therefore be considered as streambank erosion. |
The optional drainage layer can be used to represent this direct connection of the gully to the stream. To obtain the total sediment export, gully contribution (known from other studies) is added to the subwatershed level output. |
Gully erosion (disconnected from the stream) |
If the gully is disconnected, eroded soil is deposited on land and may eventually reach the stream together with rill-inter-rill sediment. |
The additional soil loss can be added in post-processing. If a map of gullies and their estimated contribution is available (e.g. Vigiak 2011, Wilkinson 2014), it may be used with the SDR layer to compute the sediment export. |
Streambank erosion |
Bank erosion is a function of soil erodibility and shear stress, i.e. the hydraulic forces applied to the channel. Detailed knowledge the bank height, and soil resistance to friction are necessary to inform a simple process-based model of these process. |
The additional soil loss can be added in post-processing. If local information is available on the magnitude of streambank erosion (e.g. per unit length of stream), users can compute this values for each subwatershed based on the stream layer generated by the model, and then add the value to the InVEST sediment export values. |
Landslides |
Landslides occur in areas of high rainfall intensity and unstable slopes. |
The additional soil loss can be analyzed using landslide-specific models and methods outside of InVEST. |
Instream deposition |
Larger particles tend to be deposited in the stream before reaching the catchment outlet |
A simple method based on the Rouse number is proposed by Pelletier (2012), which relates deposition probability to the channel slope and settling velocity. |
If you are interested in modeling in-stream processes of sediment deposition or erosion, two possibilities are CASCADE (Schmitt 2016) or Czuba 2018. Both modeling frameworks are open source, and are good if you are interested in entire river networks. If you are more interested in deposition/erosion for a smaller channel section, one option is BASEMENT (https://basement.ethz.ch/).
References¶
Bhattarai, R., Dutta, D., 2006. Estimation of Soil Erosion and Sediment Yield Using GIS at Catchment Scale. Water Resour. Manag. 21, 1635–1647.
Borselli, L., Cassi, P., Torri, D., 2008. Prolegomena to sediment and flow connectivity in the landscape: A GIS and field numerical assessment. Catena 75, 268–277.
Cavalli, M., Trevisani, S., Comiti, F., Marchi, L., 2013. Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments. Geomorphology 188, 31–41.
Czuba, J.A., 2018. A Lagrangian framework for exploring complexities of mixed-size sediment transport in gravel-bedded river networks. Geomorphology 321, 146–152. https://doi.org/10.1016/j.geomorph.2018.08.031
Desmet, P.J.J., Govers, G., 1996. A GIs procedure for automatically calculating the USLE LS factor on topographically complex landscape units. J. Soi 51, 427–433.
De Vente J, Poesen J, Verstraeten G, Govers G, Vanmaercke M, Van Rompaey, A., Boix-Fayos C., 2013. Predicting soil erosion and sediment yield at regional scales: Where do we stand? Earth-Science Rev. 127 16–29
FAO, 2006. Guidelines for soil description - Fourth edition. Rome, Italy.
Robert Griffin, Adrian Vogl, Stacie Wolny, Stefanie Covino, Eivy Monroy, Heidi Ricci, Richard Sharp, Courtney Schmidt, Emi Uchida, 2020. “Including Additional Pollutants into an Integrated Assessment Model for Estimating Nonmarket Benefits from Water Quality,” Land Economics, University of Wisconsin Press, vol. 96(4), pages 457-477. DOI: 10.3368/wple.96.4.457
Hamel, P., Chaplin-Kramer, R., Sim, S., Mueller, C. 2015. A new approach to modeling the sediment retention service (InVEST 3.0): Case study of the Cape Fear catchment, North Carolina, USA. Science of the Total Environment 524–525 (2015) 166–177.
Hughes, A.O., Prosser, I.P., 2003. Gully and Riverbank erosion mapping for the Murray-Darling Basin. Canberra, ACT.
Jamshidi, R., Dragovich, D., Webb, A.A., 2013. Distributed empirical algorithms to estimate catchment scale sediment connectivity and yield in a subtropical region. Hydrol. Process.
Lopez-vicente, M., Poesen, J., Navas, A., Gaspar, L., 2013. Predicting runoff and sediment connectivity and soil erosion by water for different land use scenarios in the Spanish Pre-Pyrenees. Catena 102, 62–73.
Oliveira, A.H., Silva, M.A. da, Silva, M.L.N., Curi, N., Neto, G.K., Freitas, D.A.F. de, 2013. Development of Topographic Factor Modeling for Application in Soil Erosion Models, in: Intechopen (Ed.), Soil Processes and Current Trends in Quality Assessment. p. 28.
Pelletier, J.D., 2012. A spatially distributed model for the long-term suspended sediment discharge and delivery ratio of drainage basins 117, 1–15.
Renard, K., Foster, G., Weesies, G., McCool, D., Yoder, D., 1997. Predicting Soil Erosion by Water: A Guide to Conservation Planning with the revised soil loss equation.
Renard, K., Freimund, J., 1994. Using monthly precipitation data to estimate the R-factor in the revised USLE. J. Hydrol. 157, 287–306. Roose, 1996. Land husbandry - Components and strategy. Soils bulletin 70. Rome, Italy.
Schmitt, R.J.P., Bizzi, S., Castelletti, A., 2016. Tracking multiple sediment cascades at the river network scale identifies controls and emerging patterns of sediment connectivity. Water Resour. Res. 3941–3965. https://doi.org/10.1002/2015WR018097
Sougnez, N., Wesemael, B. Van, Vanacker, V., 2011. Low erosion rates measured for steep , sparsely vegetated catchments in southeast Spain. Catena 84, 1–11.
Vigiak, O., Borselli, L., Newham, L.T.H., Mcinnes, J., Roberts, A.M., 2012. Comparison of conceptual landscape metrics to define hillslope-scale sediment delivery ratio. Geomorphology 138, 74–88.
Wilkinson, S.N., Dougall, C., Kinsey-Henderson, A.E., Searle, R.D., Ellis, R.J., Bartley, R., 2014. Development of a time-stepping sediment budget model for assessing land use impacts in large river basins. Sci. Total Environ. 468-469, 1210–24.