InVEST +VERSION+ documentation

Seasonal water yield model

Summary

There is a high demand for tools estimating the effect of landscape management on water supply service (e.g. for irrigation, domestic use, hydropower production). While the InVEST annual water yield model provides an estimate of total water yield for a catchment, many applications require knowledge of seasonal flows, especially during the dry season. This requires the understanding of hydrological processes in a catchment, in particular the partitioning between quick flow (occurring during or shortly after rain events) and baseflow (occurring during dry weather). In highly seasonal climates, baseflow is likely to provide greater value than the quick flow, unless significant storage (e.g., a large reservoir) is available. The InVEST seasonal water yield model seeks to provide guidance regarding the contribution of land parcels to the generation of both baseflow and quick flow. The model computes spatial indices that quantifies the relative contribution of a parcel of land to the generation of both baseflow and quick flow. Currently, there are no quantitative estimates of baseflow (only the relative contributions of pixels); a separate tool is in development to address this question.

Introduction: value attribution for individual parcels in a landscape

Understanding the effect of landscape management on seasonal flow is of critical importance for watershed management. The contribution of a given parcel to streamflow depends on a number of environmental factors including climate, soil, vegetation, slope, and position along the flow path (determining if the pixel may receive water from upslope or if water recharged may later be evapotranspired).

Water flowing across the landscape is either evaporated, transpired, withdrawn by a well, or leaves the watershed as deep groundwater flow or streamflow. If we consider an individual pixel, and its value with respect to water yield, we can consider two approaches:

  • The first gives credit to the net amount of water generated on a pixel as equal to the incoming precipitation minus the losses to evapotranspiration on that pixel. In this scheme, it is possible for actual evapotranspiration to be greater than precipitation if water is supplied to the site from upgradient. Thus, the net generation could be negative. This approach pays no heed to the eventual disposition of that water generated on that pixel; that is, it does not consider whether the water actually shows up as streamflow or is evaporated or withdrawn somewhere along its path.
  • The second approach gives credit to the water from a parcel that actually shows up as streamflow. Thus, if a parcel generates water that is later evaporated, the contribution is considered to be nil.

The former approach puts greater emphasis on the land-use and land-cover of a site, since the focus is on net generation from that pixel. It accounts for the subsidy of water from upslope pixels, but does not consider downgradient effects. It represents a potential to generate streamflow (not an actual generation of flow).

The second approach puts more emphasis on the topographic position of a pixel, as that will determine the potential for water generated on that pixel to be consumed before becoming streamflow. It represents the actual streamflow generated by a pixel. Since actual streamflow cannot be less than zero, this approach, unlike the first, will result in indices that are greater than or equal to zero.

We use these concepts to develop a set of three indices, one for quickflow, one for recharge (which represents the ‘potential baseflow’), and one for actual baseflow. Here, baseflow is defined as the generation of streamflow with watershed residence times of months to years, while quick flow represents the generation of streamflow with watershed residence times of hours to days.

How it works

QF is calculated with a CN-based approach, where we use the mean event depth, \(\frac{P_{i,m}}{n_{i,m}}\) and assume an exponential distribution of daily precipitation depths on days with rain,

\[f\left( p \right) = \frac{1}{a_{i,m}}exp\left( - \frac{p}{a_{i,m}} \right)\]

Where \(a_{i,m} = \frac{P_{i,m}}{n_{m}}/25.4\) and

  • \(a_{i,m}\) is the mean rain depth on a rainy day at pixel i on month m [in],
  • \(n_{i,m}\) is the number of events at pixel i in month m [-],
  • \(P_{i,m}\) is the monthly precipitation for pixel i at month m [mm].

In the boundary case, a stream cell’s QF is set to the precipitation at that cell

\[\text{QF}_{stream,m} = \ P_{stream,m}\]

otherwise it can be shown from the exponential distribution that the monthly runoff \(\text{QF}_{i,m}\) is

\[\text{QF}_{i,m} = n_{m} \times \left( \left( a_{i,m} - S_{i} \right)\exp\left( - \frac{0.2S_{i}}{a_{i,m}} \right) + \frac{S_{i}^{2}}{a_{i,m}}\exp\left( \frac{0.8S_{i}}{a_{i,m}} \right)E_{1}\left( \frac{S_{i}}{a_{i,m}} \right) \right) \times \left( 25.4\ \left\lbrack \frac{\text{mm}}{\text{in}} \right\rbrack \right)\]

[1]

where

  • \(S_{i} = \frac{1000}{\text{CN}_{i}} - 10\) [in]
  • \(\text{CN}_{i}\) is the curve number for pixel i [in:sup:-1], tabulated, a function of the local LULC, and soil type (see Appendix I for a template of this table),
  • and \(E_{1}\) is the exponential integral function, \(E_{1}(t) = \int_{1}^{\infty}{\frac{e^{- t}}{t}\text{dt}}\).

Thus the annual quick flow \(\text{QF}_{i}\), can be calculated from the sum of monthly \(\text{QF}_{i,m}\) values,

\[\text{QF}_{i} = \sum_{m = 1}^{12}{QF_{i,m}}\]

[2]

Local recharge (L)

The local recharge, or potential contribution to baseflow, of a pixel is computed from the local water balance. Local recharge can be negative if a pixel uses upgradient subsidies to satisfy the energy demand. The local recharge index is computed on an annual time scale, but uses values derived from monthly water budgets.

For a pixel i, the local recharge derived from the annual water budget is (Figure 1):

\[L_{i} = P_{i} - \text{QF}_{i} - \text{AET}_{i}\]

[3]

Where annual actual evapotranspiration AET is the sum of monthly AET:

\[\text{AET}_{i} = \sum_{\text{months}}^{}\text{AET}_{i,m}\]

[4]

For each month, \(\text{AET}_{i,m}\) is either limited by the demand (PET) or by the available water:

\[\text{AET}_{i,m} = min(\text{PET}_{i,m}\ ;\ P_{i,m} - \text{QF}_{i,m} + \alpha_{m}\beta_{i}L_{sum.avail,i})\]

[5]

Where \(\text{PET}_{i,m}\) is the monthly potential evapotranspiration,

\[\text{PET}_{i,m} = K_{c,i,m} \times ET_{0,i,m}\]

[6]

\(L_{sum.avail,i}\) is recursively defined by (Figure 2),

\[L_{sum.avail,i} = \sum_{j \in \{ neighbor\ pixels\ draining\ to\ pixel\ i\}}^{}{p_{\text{ij}} \cdot \left( L_{avail,j} + L_{sum.avail,j} \right)}\]

[7]

where \(p_{\text{ij}}\ \in \lbrack 0,1\rbrack\) is the proportion of flow from cell i to j, and \(L_{avail,i}\) is the available recharge to a pixel, which is 0 whenever \(L_{i}\) is negative, and a proportion \(\gamma\) of \(L_{i}\) when it is positive (see below for definition of \(\gamma\)):

\[L_{avail,i}\ = min(\gamma L_{i},L_{i})\]

[8]

In the above:

  • \(P_{i}\) and \(P_{i,m}\) are the annual and monthly precipitation, respectively [mm]
  • \(\text{QF}_{i}\) and \(\text{QF}_{i,m}\) are the quickflow indices, defined above [mm]
  • \(ET_{0,i,m}\) is the reference evapotranspiration for month m [mm]
  • \(K_{c,i,m}\) is the monthly crop factor for the pixel’s LULC [-]
  • \(\alpha_{m}\) is the fraction of upslope annual available recharge that is available in month m (default is 1/12)
  • \(\beta_{i}\) is the fraction of the upgradient subsidy that is available for downgradient evapotranspiration (default is 1; see Appendix II for more insights)
  • γ is the fraction of pixel recharge that is available to downgradient pixels (default is 1)

Attribution of recharge

The total baseflow, Qb (in mm), is the average of the contributing local recharges (negative or positive) in the catchment,

\[Q_{b} = \frac{\sum_{k \in \left\{ \text{pixels in catchment} \right\}}^{}L_{k}}{n_{\text{pixels in catchment}}}\]

[9]

Attribution value to a pixel is the relative contribution of L to the baseflow:

\[V_{R,i} = \frac{L_{i}}{{Q_{b} \times n}_{\text{pixels in catchment}}}\]

[10]

_images/fig1.png

Figure 1. Water balance at the pixel scale to compute the local recharge (Eq. 3).

_images/fig2.png

Figure 2. Routing at the hillslope scale to compute actual evapotranspiration (based on pixel’s climate variables and the upslope contribution, see Eq. 5) and baseflow (based on B:sub:`sum`, the flow actually reaching the stream, see Eq. 11-14)

The baseflow index represents the actual contribution of a pixel to baseflow (i.e. water that reaches the stream). If the local recharge is negative, then the pixel did not contribute to baseflow so B is set to zero. If the pixel contributed to groundwater recharge, then B is a function of the amount of flow leaving the pixel and of the relative contribution to recharge of this pixel.

For a parcel that is not adjacent to the stream channel, the cumulative baseflow, \(B_{sum,i}\), is proportional to the cumulative baseflow leaving the adjacent downgradient parcels minus the cumulative baseflow that was generated on that same downgradient parcel (Figure 2):

\[\begin{split}B_{sum,i} = L_{sum,i}\sum_{j \in \{\text{cells to which cell i pours}\}}^{}\begin{Bmatrix} p_{\text{ij}}\left( 1 - \frac{L_{avail,j}}{L_{sum,j}} \right)\frac{B_{sum,j}}{L_{sum,j} - L_{j}}\ \text{if }j\text{ is a nonstream pixel} \\ p_{\text{ij}}\ \text{if }j\text{ is a stream pixel} \\ \end{Bmatrix}\end{split}\]

[11]

At the watershed outlet (or at any parcel adjacent to the stream), the sum of baseflow generation \(B_{sum,i}\) over all upgradient parcels is equal to the sum of local generation over the same parcels (because there is no further opportunity for the slow flow to be consumed before reaching the stream):

\[B_{sum,outlet} = L_{sum,outlet}\]

[12]

where \(L_{sum,i}\) is the cumulative upstream recharge defined by

\[L_{sum,i} = L_{i} + \sum_{j,\ all\ pixels\ draining\ to\ pixel\ i}^{}{L_{sum,j} \cdot p_{\text{ji}}}\]

[13]

and the baseflow, \(B_{i}\) can be directly derived from the proportion of the cumulative baseflow leaving cell i, with respect to the available recharge to the upstream cumulative recharge:

\[B_{i} = max\left(B_{sum,i} \cdot \frac{L_{i}}{L_{sum,i}}, 0\right)\]

[14]

Advanced model options

One model input is the number of rain events per month, which is entered as a .csv table with one number for each month of the year. This assumes that there is one such number for the whole watershed, which may not be true for large areas or small areas with very variable precipitation fields.

To represent variability in the number of rain events, it is possible to enter a map of climate zones, and associated number of rain events for each zone.

Inputs

Name Description Type
Climate zone table Table with the number of rain events per months and climate zones. Column names: cz_id, representing climate zones numbers, integers found in the Climate zone raster, followed by 3-letter month names, i.e. “jan”,…, “dec” .csv table of integers
Climate zone Map of climate zones identified by an integer Raster of integers

The model computes sequentially the local recharge layer, and then the baseflow layer based on the first one. The local recharge layer could be obtained from a different model (e.g, RHESSys)

To compute baseflow contribution based on their own recharge layer, it is possible to bypass the first part of the model and enter directly a map of local recharge.

Inputs

Name Description Type
Local recharge Raster with the local recharge obtained from a different model (in mm) Raster of decimals

The alpha parameter represents the temporal variability in the contribution of upslope available water to evapotranspiration on a pixel. In the default parameterization, its value is set to 1/12, assuming that the soil buffers water release and that the monthly contribution is exactly one 12th of the annual contribution.

To allow upslope subsidy to be temporally variable, the user can enter the monthly αm values, in the same table as the rain events table.

Inputs

Name Description Type
Rain events table The rain events table is a model input for the default run (see below). One additional column with header “alpha” is needed to run this advanced option. .csv table

Data needs

Name Description Type
\[P_{i,m}\]
Maps of monthly precipitation (mm) Folder of 12 rasters. Rasters’ names should end with the month number (e.g. “Precip_1.tif”)
\[\text{ET}_{0,m}\]
Maps of monthly reference evapotranspiration (mm) Folder of 12 rasters. Rasters’ names should end with the month number (e.g. “ET_1.tif”)
DEM Digital elevation model Raster of decimals
LULC Map of LULC Raster of integers
Soil group Map of SCS soil hydrologic groups (A, B, C, or D), used in combination of the LULC map to compute the CN map. Values are entered as integers, with 1, 2, 3, and 4, corresponding to groups A, B, C, and D, respectively. Raster of integers
AOI/ Watershed Shapefile delineating the boundary of the area(s) of interest, or watershed to be modeled Shapefile (can be polyshape)
Biophysical table

Table comprising, for each LULC type:

  • CN for each soil type
  • Monthly Kc values
.csv file with column names: CN_A, CN_B, CN_C, CN_D, Kc_1, …, Kc_12
Rain events table Table with 12 values of rain events per month. A rain event is defined as >0.1mm (USGS: http://drought.unl.edu/MonitoringTools/USRainDaysandDryDays.aspx) .csv file with column names: month and events
Threshold flow accumulation 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 in the DEM. Integer
\(\alpha_{m}\), \(\beta_{i}\), γ

Model parameters used for research purposes. Default values are:

\(\alpha_{m} = 1/12\), \(\beta_{i} = 1\), γ=1

Decimal

Data sources and guidance for parameter selection

Name Source
\[P_{i,m}\]

Global monthly precipitation data can be obtained from the WorldClim dataset: http://www.worldclim.org/

Alternatively, rasters can be extrapolated from rain gauges with monthly data.

\[\text{ET}_{0,m}\]
Global monthly reference evapotranspiration may be obtained from the CGIAR CSI dataset (based on WorldClim data): http://www.cgiar-csi.org/data/global-aridity-and-pet-database
DEM

DEM data is available for any area of the world, although at varying resolutions.

Free raw global DEM data is available from:

Alternatively, it may be purchased relatively inexpensively at sites such as MapMart (www.mapmart.com).

The DEM resolution may be a very important parameter depending on the project’s goals. For example, if decision makers need information about impacts of roads on ecosystem services then fine resolution is needed. The hydrological aspects of the DEM used in the model must be correct. Because the model requires that all pixels have a flow direction (according to the D-infinity flow algorithm (Tarboton, 1997)), the DEM may need to be filled to remove sinks. Multiple passes of the ArcGis Fill tool, or Qgis Wang&Liu Fill algorithm (SAGA library) have shown good results.

LULC

A key component for all water models is a spatially continuous landuse / land cover raster grid. That is, within a watershed, all landuse / land cover categories should be defined. Gaps in data will create errors. Unknown data gaps should be approximated. Global land use data is available from:

Data for the U.S. for 1992 and 2001 is provided by the EPA in their National Land Cover Data product: http://www.epa.gov/mrlc/.

The simplest categorization of LULCs on the landscape involves delineation by land cover only (e.g., cropland, temperate conifer forest, prairie). Several global and regional land cover classifications are available (e.g., Anderson et al. 1976), and often detailed land cover classification has been done for the landscape of interest.

A slightly more sophisticated LULC classification involves breaking relevant LULC types into more meaningful types. For example, agricultural land classes could be broken up into different crop types or forest could be broken up into specific species. The categorization of land use types depends on the model and how much data is available for each of the land types. Users should only break up a land use type if it will provide more accuracy in modeling. For instance, for the sediment model the user should only break up ‘crops’ into different crop types if they have information on the difference in soil characteristics between crop management values.

Soil group

Soil groups are determined from hydraulic conductivity and soil depths.

FutureWater has created a global map of hydraulic conductivity available at: http://www.futurewater.eu/2015/07/soil-hydraulic-properties/

To convert hydraulic conductivity to soil hydrologic group, Table 1 below can be used.

Otherwise, one can find guidance online, e.g.: www.bwsr.state.mn.us/outreach/eLINK/Guidance/HSG_guidance.pdf

AOI/ Watershed

To delineate watersheds, users can use the InVEST tool DelineateIT

Alternatively, a number of watershed maps are available online, e.g. HydroBASINS: http://hydrosheds.org/

Biophysical table
  • CN can be obtained from the USDA handbook: (NRCS-USDA, 2007 Chap. 9)
  • Monthly Kc values can be obtained from the FAO guidelines: (Allen et al., 1998)

For water bodies and wetlands that are connected to the stream, CN can be set to 99 (i.e. assuming that those pixels rapidly convey quickflow)

Note: when the focus is on potential flood effects, CN may be selected to reflect wet antecedent runoff conditions: CN values should then be converted to ARC-III conditions, as per Chap 10 in NRCA-USDA guidelines (2007)

Rain events table

The average number of monthly rain events can be obtained from local climate statistics (Bureau of Meteorology) or other online resources (eg.http://www.yr.no/, http://wcatlas.iwmi.org). The World Bank also provides maps with precipitation statistics: http://data.worldbank.org/developers/climate-data-api

Climate zones from: http://koeppen-geiger.vu-wien.ac.at/present.htm

(to delineate a reasonable number of zones)

Threshold flow accumulation

Needs to be adjusted based on local stream maps.

Rule of thumb: contribution area of 1km2 (threshold needs to be calculated based on pixel area)

Stream maps can be obtained from HydroSHEDS:

\[\alpha_{m}\]
Default=1/12. See Appendix
\[\beta_{i}\]
Default=1. See Appendix
γ Default =1. See Appendix

Table 1: Criteria for assignment of hydrologic soil groups (NRCS-USDA, 2007 Chap. 7)

  Group A Group B Group C Group D
Saturated hydraulic conductivity of the least transmissive layer when a water impermeable layer exists at a depth between 50 and 100 centimeters >40 μm/s [40;10] μm/s [10;1] μm/s <1 μm/s (or depth to impermeable layer<50cm or water table<60cm)
Saturated hydraulic conductivity of the least transmissive layer when any water impermeable layer exists at a depth greater than 100 centimeters >10 μm/s [4;10] μm/s [0.4;4] μm/s <0.4 μm/s

Interpreting outputs

  • CN (raster): Map of CN values
  • QF (raster): Map of quickflow QF values [mm]
  • L (raster): Map of local recharge \(L\) values [mm]
  • L_avail (raster): Map of available local recharge \(L_{\text{avail}}\) , i.e. only positive L values [mm]
  • B (raster): Map of baseflow \(B\) values [mm], the contribution of a pixel to slow release flow (which is not evapotranspired before it reaches the stream)
  • B_sum (raster): Map of \(B_{\text{sum}}\)values [mm], the flow through a pixel, contributed by all upslope pixels, that is not evapotranspirated before it reaches the stream
  • L_sum (raster): Map of \(L_{\text{sum}}\) values [mm], the flow through a pixel, contributed by all upslope pixels, that is available for evapotranspiration to downslope pixels
  • L_sum_avail (raster): Map of \(L_{\text{sum.avail}}\) values [mm], the available water to a pixel, contributed by all upslope pixels, that is available for evapotranspiration by this pixel
  • Q_b (decimal): Annual average baseflow [mm]
  • V_Ri (raster): Map of the values of recharge (contribution, positive or negative, to the total recharge

References:

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration - Guidelines for computing crop water requirements, FAO Irrigation and drainage paper 56. Rome, Italy.

NRCS-USDA, 2007. National Engineering Handbook. United States Department of Agriculture, http://www.nrcs.usda.gov/wps/portal/nrcs/detailfull/national/water/?cid=stelprdb1043063.

Appendix: \({\mathbf{\alpha},\mathbf{\beta}}_{\mathbf{i}},\)and γ parameters definition and alternative values

\(\alpha\) and \(\beta_{i}\) represent the fraction of annual recharge from upgradient parcels that is available to a downgradient pixel for evapotranspiration in a given month. The product \(\alpha \times \beta_{i}\) is expected to be <1 since some water from upslope may be unavailable, either when it follows deep flowpaths or when the timing of supply and (evapotranspirative) demand is not right.

\(\alpha\) is a function of precipitation seasonality: recharge from a given month can be used by downslope areas during later months, depending on the subsurface travel times. In the default parameterization, its value is set to 1/12, assuming that the soil buffers water release and that the monthly contribution is exactly one 12th of the annual contribution. An alternative assumption is to set values to the antecedent monthly precipitation values, relative to the total precipitation: Pm-1/P:sub:annual

\(\beta_{i}\) is a function of local topography and soils: for a given amount of upslope recharge, the amount of water used by a pixel is a function of the storage capacity. It also depends on the characteristics of the upslope area: the use of the upgradient subsidy is conditioned by the shape and area of the contribution area (i.e. the recharge from the pixel just above the pixel of interest is less likely to be lost than the pixels much further away)

In the default parameterization, \(\beta\) is set to 1 for all pixels. One alternative is be to set \(\beta_{i}\) as TI, the topographic wetness index for a pixel, defined as \(ln(\frac{A}{\text{tan}\beta}\)) (or other formulation including soil type and depth).

γ represents the fraction of pixel recharge that is available to downgradient pixels. It is a function of soil properties and possibly topography (e.g. with very permeable soils, the value of . In the default parameterization, γ is constant over the landscape and plays a role similar to \(\alpha\).

In practice

The options above are provided mainly for research purposes. In practice, we suggest that for highly seasonal climates, alpha should be set to the antecedent monthly precipitation values, relative to the total precipitation: Pm-1/P:sub:annual

Then, we offer two options to address the uncertainty around the parameter values:

  1. Verification of actual evapotranspiration with observations

The model outputs the actual evapotranspiration at the annual time scale: users can adjust parameters to meet observed actual evapotranspiration (e.g. from MODIS, http://www.ntsg.umt.edu/project/mod16).

  • If AET_mod>AET_obs, the model overpredicts evapotranspiration, which can be corrected by: reducing Kc values, or reducing gamma values, and/or beta values (so less water is available for each pixel).
  • If AET_mod<AET_obs, the model underpredicts evapotranspiration, which can be corrected by: increasing Kc values (and increasing gamma or beta values if they are not at their maximum of 1).

If monthly values of AET are available, a finer calibration can be performed by changing the seasonal parameter alpha.

  1. Ensemble modeling

The model can be run under different assumptions and the outputs compared to estimate the effect of parameter error. Parameter ranges can be determined from assumptions about the proportion of upslope subsidy available to a given pixel; they can be set to the maximum bounds (0 and 1) for preliminary results.