Chapter 19 Watershed hydrology

This chapter describes the distributed watershed water balance sub-model implemented in functions spwb_land(), growth_land() and fordyn_land() of package medfateland (Fig. 1.3).

19.1 Design principles

Watersheds are described in raster (i.e. gridded) mode, each cell representing a patch of vegetation (or another land cover) in a catchment. Land cover types are:

  1. wildland: forests, shrublands or grasslands.
  2. agriculture: agricultural lands.
  3. rock: rock outcrops.
  4. artificial: urban areas.
  5. water: water bodies.

Watershed water balance simulations extend the simulation of soil water balance to the landscape scale. Hence, the design of most vertical hydrological processes is the same as those of forest water balance. However, additional water compartments and processes are represented in watershed water balance simulations. Each cell in the watershed may include the following: a snow-pack compartment, one or several soil layers (including a rocky layer down to several meters depth), and a groundwater compartment assumed to be beyond the reach of plant roots. Similar to other models such as TETIS (Francés et al. 2007), three lateral flows are considered between adjacent cells:

  1. Overland surface flows from upper elevation cells.
  2. Lateral saturated soil flows (i.e. interflow) between adjacent cells.
  3. Lateral groundwater flow (i.e. baseflow) between adjacent cells.

Overland surface flows are modeled following T-HYDRO Ostendorf & Reynolds (1993). Subsurface flows, including both soil saturated flow and groundwater flow, are modeled following the kinematic wave approach of DSHVM Wigmosta et al. (1994).

Split-parameter parametrization following Francés et al. (2007).

19.2 State variables

Distributed simulation of rainfall-runoff processes implies that state variables are defined for each cell of the watershed. For this kind of simulations, state variables describing the water content of different compartments are the most important. In the case of wildland and agriculture cells these are:

  • \(W_s\), the proportion of soil moisture in relation to field capacity for each soil layer \(s\) in the cell.
  • \(S_{snow}\) the snow water equivalent (mm) of the snow pack storage over the cell surface.
  • \(S_{aquifer}\) the water content (mm) in the cell’s aquifer beyond the reach of plant roots.

Additional cell state variables in wildland cells concern the water status (or other state variables) of plant cohorts in the cell, and they were described in chapters 3 and 7. In the case of cells of other land cover types (i.e. rock, artificial or water) soils are not defined, and neither are soil state variables.

19.3 Water balance

In distributed watershed simulations water balance can be defined both at the cell level and at the watershed level.

19.3.1 Cell level water balance

At the cell level, daily variations in water content can occur in the snowpack, soil or groundwater compartments:

  1. Variations in snowpack water equivalent content follows eq. (3.2).
  2. Variations in soil water content (\(\Delta{V_{soil}}\)) need to account for additional flows and are summarized as (compare to (3.1)): \[\begin{equation} \Delta{V_{soil}} = Pr + Sm + Ro + Ad - In - Ru - Dd - Es - Tr + \Delta{S_{sat}} \tag{19.1} \end{equation}\] where \(Pr\) is precipitation as rainfall, \(Sm\) is snowmelt, \(Ro\) is surface runon water entering the cell from neighboring cells at higher elevation, \(Ad\) is aquifer discharge, \(In\) is rainfall interception loss, \(Ru\) is cell surface runoff, \(Dd\) is deep drainage towards the aquifer, \(Es\) is soil evaporation, \(Tr\) is plant transpiration and \(\Delta{S_{sat}}\) is the variation in soil water content derived from the balance between saturated soil lateral inputs and outputs.
  3. Variations in the aquifer water content of a cell are summarized by: \[\begin{equation} \Delta{S_{aquifer}} = Dd - Ad - \Delta{S_{base}} \tag{19.2} \end{equation}\] where \(\Delta{S_{base}}\) is the variation in water content of the aquifer derived from the balance between groundwater lateral inputs and outputs.

19.3.2 Watershed water balance

At the watershed level, separate water balances can again be defined for the average water content of snow pack, soil or groundwater compartments. These result from averaging water balances across cells. Additionally, a water balance is defined regarding the overall water content in the watershed.

  1. Changes in the average snow pack water equivalent over cells is the result of balancing precipitation as snow (\(Ps\)) and snow melt (\(Sm\)), both flows averaged over cells: \[\begin{equation} \Delta{\hat{S}_{snow}} = \hat{Ps} - \hat{Sm} \end{equation}\] where \(\hat{Ps}\) and \(\hat{Sm}\) are the average snow fall and snow melt over cells.
  2. Changes in the average soil moisture (\(\Delta{\hat{V}_{soil}}\)) are the result of pooling soil inputs and outputs over cells, which yields: \[\begin{equation} \Delta{\hat{V}_{soil}} = \hat{Pr} + \hat{Sm} + \hat{Ad} + \hat{Ro} - \hat{Ru} - \hat{In} - \hat{Dd} - \hat{Es} - \hat{Tr} \end{equation}\] where \(\hat{Pr}\), \(\hat{Ad}\), \(\hat{Ro}\), \(\hat{Ru}\), \(\hat{In}\), \(\hat{Dd}\), \(\hat{Es}\) and \(\hat{Tr}\) are the average (over cells) of precipitation (including rain and snow), aquifer discharge, lateral surface water input (runon), runoff, rainfall interception loss, deep drainage, soil evaporation and plant transpiration, respectively. Lateral saturated soil flows are not included as they cancel out at the watershed level.
  3. Changes in the average of cell aquifer water content (\(\Delta{\hat{S}_{aquifer}}\)) are the result of balancing deep drainage (\(Dd\)) from soils and aquifer discharge (\(Ad\)): \[\begin{equation} \Delta{\hat{S}_{aquifer}} = \hat{Dd} - \hat{Ad} \end{equation}\]
  4. If we integrate the three water compartments, water balance at the watershed level is given by: \[\begin{equation} \Delta{S_{watershed}} = \hat{Pr} - \hat{In} - \hat{Es} - \hat{Tr} - \hat{Ex} \end{equation}\] where \(\Delta{S_{watershed}} = \Delta{\hat{S}_{snow}}+\Delta{\hat{V}_{soil}}+\Delta{\hat{S}_{aquifer}}\) is the change in water content in the watershed and \(\hat{Ex}\) is the water exported as runoff from cells without neighbors (i.e. catchment outlet cells).

19.4 Process scheduling

For every day to be simulated, the model performs the following steps:

  1. Calculation of soil and groundwater hydraulic heads for each cell (see 19.6.1).
  2. Calculation of soil and groundwater lateral flows between adjacent cells, according to hydraulic gradients (see 19.6.1).
  3. Apply changes in cell soil moisture content due to soil lateral flows (i.e. determine \(\Delta{S_{sub}}\) for each cell), including the possibility of a return flow to the surface if the soil becomes saturated (i.e. saturation excess flow).
  4. Apply changes in cell aquifer water content, including the possibility of discharge from the aquifer to the soil and a saturation excess flow (eq. (19.2)).
  5. Determine remaining cell flows by processing watershed cells in order of decreasing elevation (i.e. cells at higher elevation are processed before cells at lower elevation). For each cell to be processed:
    1. For wildland or agriculture cells, determine snow-pack dynamics, rainfall interception loss, infiltration and runoff, transpiration, soil evaporation processes as described in 3.4, while including saturation excess flow as well as potential runon (\(Ro\)) from upslope cells as additional water inputs. For rock and artificial cells, snow pack dynamics is still processed in these cells, but all runon or rainfall input becomes runoff to be passed downhill, and when the level of the aquifer reaches the surface any discharge also becomes runoff. In the case of water cells, snow pack dynamics is also processed, but evaporation from the water surface is not modelled. All liquid water inputs are poured onto the aquifer water content, and runoff can still occur when the aquifer reaches the surface.
    2. Distribute surface runoff among cells downhill (19.6.2). If a cell does not have downhill neighbors (i.e outlet cell) its runoff becomes water exported from the watershed.

19.5 Inputs and outputs

19.5.1 Gridded inputs

Hydrological distributed models demand large amounts of data, information and parameters in order to accurately represent the spatial variability of the main hydrological processes and weather inputs. As mentioned in design principles, watersheds are described in raster (i.e. gridded) mode, each cell representing a patch of vegetation (or another land cover) in a catchment. Rasters do not need to be (should not be) complete, since the model is intended to be used on individual catchments with defined boundaries.

Given a set of grid cells representing the watershed at a given spatial resolution, input data for spwb_land() and growth_land() should be specified using sf objects. These were described in section 2.1.2.

19.5.2 Watershed meteorology input

Simulation of watershed hydrology requires specifying weather variables for all cells and the simulation period considered. Functions spwb_land() and growth_land() allow weather to be specified in three ways:

  • If a single data frame of daily weather is provided, the model assumes that all watershed cells experience the same weather.
  • A different data frame of daily weather may be provided for each cell, using column meteo of the sf object (see 2.1.2).
  • If an interpolator object of class stars is provided, the model will call function interpolate_data() for every day to be simulated (see meteoland package documentation).

19.5.3 Watershed hydraulic correction parameters

Distributed hydrological models suffer from the need to specify many parameters. Frances et al. (2007) demonstrated that by using a proper parameter structure, it is possible to obtain an excellent automatic calibration of a distributed conceptual model while maintaining the spatial variability of parameters. In particular, they showed how to define effective model parameters for those quantities normally estimated at the point scale, so that model parameters and estimates usually suffer from being uncommensurable. Effective model parameters can in this situation be defined as a correction function of the corresponding hydrological characteristic. If the correction function takes into account the model and input errors, the temporal and spatial scale effects and also the hydrological characteristics estimation error, it is reasonable to assume the correction function for each parameter will be common to all cells within the watershed. Frances et al. (2007) suggest using a correction factor for each hydraulic parameter, in a strategy called split-parameter. From the calibration point of view, a very important consequence of the split-parameter approach is that the number of variables to be adjusted is reduced dramatically.

Our model follows Frances et al. (2007) in that correction factors are defined for three hydraulic parameters with spatial variation across the watershed:

  • \(K_{drain}\) vertical saturated soil conductivity for deep drainage.
  • \(K_{sat}\) horizontal saturated soil conductivity for soil lateral flows.
  • \(K_{base}\) bedrock conductivity for groundwater lateral flows.

The effective model parameter for each of these is defined for the spatial (depending on the grid resolution) and temporal (daily) scales of the model. Each effective parameter is split into the input value, assumed to have been estimated at a given spatial and temporal scale, and the corresponding correction factor, i.e. \(R_{drain}\), \(R_{sat}\) and \(R_{base}\). The values for these correction factors should be the target of calibration procedures, assuming that the spatial heterogeneity of effective parameters is adequately covered by the soil and bedrock spatial inputs.

19.5.4 Model outputs

Distributed hydrological models can produce very detailed spatial and temporal outputs, hence consuming lots of memory resources. While the model runs at a daily scale, functions spwb_land() and growth_land() can be asked to retain spatial outputs and summaries at coarser temporal scales (i.e. weekly, monthly, etc.). The output of functions spwb_land() and growth_land() is an S3 list with the following data elements:

Element Description
sf An object of class sf containing the geometries, final state and summary of the landscape cells.
watershed_balance Data frame with averaged (over cell) components of the watershed balance corresponding to each summary.
watershed_soil_balance Averaged (over cell) components of the watershed balance, including only cells that have defined soils (i.e. wildland and agriculture cells), corresponding to each summary.
daily_runoff A matrix with daily runoff (in m3/day) at each of the outlet cells of the landscape.

19.6 Process details

19.6.1 Subsurface flows

Both soil saturated lateral flow and groundwater lateral flow follow the saturated subsurface flow model of Wigmosta et al. (1994), which is based on a kinematic wave approximation. Each grid cell can exchange water with its eight adjacent neighbors. Local hydraulic gradients are approximated by local ground surface slopes, so that each cell will generally receive water from its upslope neighbors and discharge to its downslope neighbors. The rate of saturated flow from cell \(i\) to its downgradient neighbors is equal to (Wigmosta & Lettenmaier 1999): \[\begin{eqnarray} q_{ij} &=& - T_{i}\cdot \tan(\beta_{i,j})\cdot w_j \;\; \beta_{i,j} < 0\\ q_{ij} &=& 0 \;\; \beta_{i,j} \geq 0 \end{eqnarray}\] where \(q_{i,k}\) is the flow rate from cell \(i\) in the direction towards neighbor \(j\), \(T_{i}\) is the transmissivity at cell \(i\), \(\beta_{i,j}\) is the water table slope in the direction of \(j\). The transmissivity function for soil saturated lateral flow in cell \(i\) (\(T_{sat,i}\)) is specified using: \[\begin{equation} T_{sat,i} = \frac{R_{sat}\cdot K_{sat,i} \cdot Z_{soil,i}}{n_{sat}} \cdot (1 - Z_{wt,i}/Z_{soil,i})^{n_{sat}} \end{equation}\] where \(K_{sat,i}\) is the saturated conductivity of the first (top) soil layer in cell \(i\), \(Z_{soil,i}\) is the soil depth at cell \(i\), \(Z_{wt,i}\) is the water table depth at cell \(i\) and \(n_{sat}\) is a power exponent. \(R_{sat}\) is a correction factor (a watershed parameter) for the fact that \(K_{sat,i}\) is assessed from soil samples, whereas the parameter in the model is defined at the scale of a grid cell (Francés et al. 2007). The transmissivity function for groundwater lateral flow in cell \(i\) (\(T_{base,i}\)) is analogous: \[\begin{equation} T_{base,i} = \frac{R_{base}\cdot K_{base,i} \cdot Z_{bedrock,i}}{n_{base}} \cdot (1 - Z_{aquifer,i}/Z_{bedrock,i})^{n_{base}} \end{equation}\] where \(K_{base,i}\) is the bedrock conductivity at cell \(i\), \(Z_{bedrock,i}\) is the depth of the unaltered bedrock, \(Z_{aquifer,i}\) is the depth of the aquifer water table and \(n_{base}\) is again a power exponent.

Saturated lateral flow inputs and outputs for a given cell \(i\) are balanced to determine increase or decrease in soil moisture: \[\begin{equation} \Delta{S_{sat,i}} = \sum_{j}{q_{sat, ji}} - \sum_{j}{q_{sat, ij}} \end{equation}\] and the same occurs for groundwater lateral flows and aquifer balance: \[\begin{equation} \Delta{S_{base,i}} = \sum_{j}{q_{base, ji}} - \sum_{j}{q_{base, ij}} \end{equation}\]

19.6.2 Overland flows

To simulate surface runoff routing from one cell to the other, the approach of Ostendorf & Reynolds (1993) is used, as in SIERRA (Mouillot et al. 2001). Overland water lateral transport for a given day occurs instantaneously (i.e. no velocities are calculated) and depends on topography only. The model determines cell neighbors following the queen rule (up to eight neighbors per cell). The proportion of surface water runoff of cell \(i\) that will be added to the infiltration input (runon) of a neighboring cell \(j\) is Ostendorf & Reynolds (1993): \[\begin{equation} q_{ij} = \frac{\Delta z_{ij}/L_{ij}}{\sum_{j}{\Delta z_{ij}/L_{ij}}} \end{equation}\] if \(\Delta z_{ij} = z_i - z_j > 0\), that is, if the difference in elevation between the two cells is positive (i.e. if \(z_j < z_i\)). Otherwise there is no overland transport from \(i\) to \(j\), i.e. \(q_{ij} = 0\). \(L_{ij}\) indicates the distance between cell \(i\) and \(j\) (which depends on cell size and on whether the neighboring cell \(j\) is diagonal to cell \(i\)). The summation of the denominator is done only for neighbors at lower elevation, so that \(\sum_{i}{q_{ij}} = 1\). The table of \(q_{ij}\) values is calculated when initializing distributed watershed objects.

Every day, cells are processed in order from higher to lower elevation. After the daily water balance of a given cell \(i\), water runoff \(Ru_i\) is divided among the neighboring cells at lower elevation. The runon of a neighbor \(j\), \(Ro_j\) is updated as:

\[\begin{equation} Ro_j = Ro_j + Ru_i \cdot q_{ij} \end{equation}\]

Note that a given cell \(j\) can receive overland flow from more than one neighbor.

19.6.3 Watershed runoff

A special situation arises when processing cells that do not have downhill neighbors defined (i.e. where \(q_{ij} = 0\) for all \(j\)), either because they are in flat surfaces or at the watershed boundary. In both cases, these cells should correspond to water bodies or streams connecting to the catchment outlet. Hence, these cells are called outlet cells in the model, and the runoff they generate becomes watershed runoff.


Francés, F., Vélez, J.I. & Vélez, J.J. (2007). Split-parameter structure for the automatic calibration of distributed hydrological models. Journal of Hydrology, 332, 226–240.
Mouillot, F., Rambal, S. & Lavorel, S. (2001). A generic process-based SImulator for meditERRanean landscApes (SIERRA): Design and validation exercises. Forest Ecology and Management, 147, 75–97.
Ostendorf, B. & Reynolds, J.F. (1993). Relationships between a terrain-based hydrologic model and patch-scale vegetation patterns in an arctic tundra landscape. Landscape Ecology, 8, 229–237.
Wigmosta, M.S. & Lettenmaier, D.P. (1999). A comparison of simplified methods for routing topographically driven subsurface flow. Water Resources Research, 35, 255–264.
Wigmosta, M., Vail, L. & Lettenmaier, D. (1994). A distributed hydrology-vegetation model for complex terrain. Water Resources Research, 30, 1665–1679.