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. Five land cover types are allowed:
- wildland: forests, shrublands or grasslands.
- agriculture: agricultural lands.
- rock: rock outcrops.
- artificial: urban areas.
- 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, cells classified as rock, artificial or water lack any soil and must be treated differently. 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 normally beyond the reach of plant roots for most landscape cells. Similar to other models such as TETIS (Francés et al. 2007), three lateral flows are considered between adjacent cells:
- Overland surface flows from upper elevation cells.
- Lateral saturated soil flows (i.e. interflow) between adjacent cells.
- Lateral groundwater flow (i.e. baseflow) between adjacent cells.
Routing of overland surface flows is modeled following T-HYDRO Ostendorf & Reynolds (1993). Subsurface flows, including both soil saturated flow and groundwater flow, are modeled following the kinematic wave approximation of DSHVM Wigmosta et al. (1994).
Split-parameter parametrization, following Francés et al. (2007), allows reducing the parameters to be calibrated for practical applications.
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) of the cell’s aquifer lying over the undisturbed bedrock.
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. In this case, only \(S_{snow}\) and \(S_{aquifer}\) are defined.
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:
- Variations in snowpack water equivalent content follows the balance between snowfall (\(Ps\)) and snow melt (\(Sm\)) eq. (3.1): \[\begin{equation} \Delta{S_{snow}} = Ps - Sm \tag{19.1} \end{equation}\]
- Variations in soil water content (\(\Delta{V_{soil}}\)) follow (3.3): \[\begin{equation} \Delta{V_{soil}} = (If + Lw + Cr) - (Dd + Se + Es + Tr_{herb} + Tr_{woody}) \tag{19.2} \end{equation}\] where the terms were defined in 3.3. Here, the lateral flows (i.e. \(Lw\)) are specially important, compared to spatially-uncoupled simulations. They are defined as \(Lw = \Delta S_{lat}\), the variation in water content derived from the balance between soil lateral inputs and outputs: \[\begin{equation} Lw = \Delta S_{lat} = S_{lat, in} - S_{lat, out} \end{equation}\]
- Variations in the aquifer water content of a cell are summarized by: \[\begin{equation} \Delta{S_{aquifer}} = Dd - Cr - Exf - Dal - \Delta{S_{base}} \tag{19.3} \end{equation}\] where \(Dd\) is deep drainage from the above-lying soil, \(Cr\) is capillary rise into the above-lying soil, \(Exf\) is the aquifer exfiltration (ocurring when the water table reaches the soil surface), \(Dal\) is drainage loss to a deeper aquifer (not connected to watershed outlets) and \(\Delta{S_{base}}\) is the variation in water content of the aquifer derived from the balance between groundwater lateral inputs and outputs. \[\begin{equation} \Delta S_{base} = S_{base, in} - S_{base, out} \end{equation}\]
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.
- 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.
- 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{If} + \hat{Cr}) - (\hat{Dd} + \hat{Se} + \hat{Es} + \hat{Tr}_{woody} + \hat{Tr}_{herb}) \end{equation}\] where \(\hat{If}\), \(\hat{Cr}\), \(\hat{Dd}\), \(\hat{Se}\), \(\hat{Es}\), \(\hat{Tr}_{woody}\) and \(\hat{Tr}_{herb}\) are the average (over cells) of infiltration, capillary rise, soil evaporation, plant transpiration and herb transpiration, respectively. Lateral saturated soil flows are not included as they cancel out at the watershed level.
- Changes in the average of cell aquifer water content (\(\Delta{\hat{S}_{aquifer}}\)) are the result of balancing averages for deep drainage (\(\hat{Dd}\)) from soils, capillarity rise to soils (\(\hat{Cr}\)), aquifer exfiltration (\(\hat{Exf}\)) and loss to a deeper aquifer (\(\hat{Dal}\)): \[\begin{equation} \Delta{\hat{S}_{aquifer}} = \hat{Dd} - \hat{Cr} - \hat{Exf} - \hat{Dal} \end{equation}\]
- 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}_{woody} + \hat{Tr}_{herb}) - \hat{Exp} - \hat{Dal} \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{Exp}\) 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:
- Calculation of soil lateral flows between adjacent cells, according to sub-steps if required by the user (see 19.6.1).
- Determine the balance of soil lateral flows (i.e. determine \(\Delta{S_{lat}}\)) for each cell.
- Calculation of groundwater lateral flows between adjacent cells (see 19.6.1).
- Apply changes in cell aquifer water content due to the balance of groundwater lateral flows (\(\Delta{S_{base}}\) in eq. (19.3)).
- Process vertical processes in cells that lack soil (i.e. rock, artificial and water). For rock and artificial cells, snow pack dynamics is still processed in these cells. In artificial cells all rainfall input or snow melt becomes runoff to be passed downhill. In rock cells, water can infiltrate directly to the aquifer, where a maximum rate is confronted with the rainfall/snowmelt input. In the case of water cells, snow pack dynamics are still processed, but all liquid water (from rain or snow melt) is poured (as deep drainage) down to the aquifer below.
- Distribute overland flows from rock or artificial cells donwhill to wildland or agriculture cells, which will receive some runon.
- Determine vertical flows in watershed cells having soil (i.e. wildland or agriculture). Cells are processed sequentially or in parallel, while accounting for subsurface flows, runon from upslope rock or artificial cells, and changes in the groundwater level, which complete eqs. (19.1) and (19.2). All processes mentioned in 3.4 are taken into account.
- Distribute overland flows from cells having soil 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.6.3).
- Apply vertical flows between soils and the aquifer below, i.e. deep drainage and capillarity rise, to the water balance of the aquifer.
- Apply drainage losses from the aquifer towards a deeper aquifer not connected to watershed outlets (19.6.4), completing eq. (19.3).
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. Each row in the sf object has to correspond to a single raster cell, but there can be raster cells that are not described in the sf object and, hence, correspond to empty cells.
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
meteoof thesfobject (see 2.1.2). - If an interpolator object of class
starsis provided, the model will call functioninterpolate_data()for every day to be simulated (see meteoland package documentation).
19.5.3 Watershed hydraulic control 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_{sat, vert}\) saturated soil conductivity for vertical fluxes into, within and out of the soil.
- \(K_{sat, lat}\) saturated soil conductivity for soil lateral (horizontal) fluxes.
- \(K_{sat, 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. Parameters \(K_{sat, vert}\) and \(K_{sat, lat}\) are split into the input soil saturated conductivity value \(K_{sat}\) (in the case of \(K_{sat, lat}\), that corresponding to the topsoil), assumed to have been estimated at a given spatial and temporal scale, and the corresponding correction factor, i.e. \(R_{vert}\) or \(R_{lat}\) . Analogously, \(K_{sat, base}\) is split between an input bedrock conductivity value (\(K_{bedrock}\), supplied as bedrock_conductivity layer) and the correction factor \(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.
The following table describes the full set of watershed hydraulic control parameters:
| Symbol | R param | Description |
|---|---|---|
| \(R_{vert}\) | R_localflow |
Multiplier of soil vertical conductivity. |
| \(R_{lat}\) | R_interflow |
Multiplier of soil lateral conductivity (eq. (19.4)). |
| \(R_{base}\) | R_baseflow |
Multiplier of groundwater lateral conductivity (eq. (19.5)). |
| \(n_{lat}\) | n_interflow |
Power exponent regulating the effect of soil saturation level on soil lateral transmissivity (eq. (19.4)). |
| \(n_{base}\) | n_baseflow |
Power exponent regulating the effect of aquifer height on groundwater lateral transmissivity (eq. (19.5)). |
num_daily_substeps |
Number of daily sub-steps for interflow calculations. | |
rock_max_infiltration |
Maximum infiltration rate for rocky cells (mm/day). | |
deep_aquifer_loss |
Maximum daily loss rate from watershed aquifer towards a deeper aquifer not connected to outlets (mm·day-1). |
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 |
|---|---|
watershed_control |
A list with the watershed control parameters used in the simulation. |
sf |
An object of class sf containing the geometries, final state and temporal summary for each of the landscape cells (values in mm/day, mm/week, mm/month or mm/year depending on summary frequency). Additionally, it can contain full water balance simulation results for specific cells targetted by the user. |
watershed_balance |
Data frame with spatially-averaged (over all watershed cells) components of the daily water balance. |
watershed_soil_balance |
Data frame with spatially-averaged components of the daily water balance, including only cells that have defined soils (i.e. wildland and agriculture cells). |
outlet_export_m3 |
A matrix with daily runoff (in m3/day) at each of the outlet cells of the landscape. |
19.6 Process details
19.6.1 Soil and groundwater lateral fluxes
Both soil saturated and groundwater lateral fluxes 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 flux rate (in \(m3/day\)) from cell \(i\) in the direction towards neighbor \(j\), \(T_{i}\) is the transmissivity (in \(m2/day\)) at cell \(i\), \(\beta_{i,j}\) is the water table slope in the direction of \(j\) and \(w_j\) is the width (in \(m\)) of cell \(j\). The transmissivity function for soil saturated lateral flux in cell \(i\) (\(T_{lat,i}\)) is specified using: \[\begin{equation} T_{sat,i} = \frac{R_{lat}\cdot K_{sat,i} \cdot Z_{soil,i}}{n_{lat}} \cdot (1 - Z_{sat,i}/Z_{soil,i})^{n_{lat}} \tag{19.4} \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_{sat,i}\) is the depth of saturated water at cell \(i\) (see 2.3.4) and \(n_{lat}\) is a power exponent regulating the responsiveness of transmissivity to the level of saturation. \(R_{lat}\) is the 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), as explained above.
The transmissivity function for groundwater lateral flux in cell \(i\) (\(T_{base,i}\)) is analogous: \[\begin{equation} T_{base,i} = \frac{R_{base}\cdot K_{bedrock,i} \cdot Z_{bedrock,i}}{n_{base}} \cdot (1 - Z_{aquifer,i}/Z_{bedrock,i})^{n_{base}} \tag{19.5} \end{equation}\] where \(K_{bedrock,i}\) is the bedrock conductivity at cell \(i\), \(R_{base}\) is the correction factor for base flows, \(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 regulating the responsiveness of transmissivity to the aquifer level.
Soil saturated lateral flux inputs and outputs for a given cell \(i\) are balanced to determine increase or decrease in soil moisture (to be incorporated as source/sink when processing vertical soil flows): \[\begin{equation} \Delta{S_{lat,i}} = S_{lat, in} - S_{lat, out} = \sum_{j}{q_{lat, ji}} - \sum_{j}{q_{lat, ij}} \end{equation}\] and the same occurs for groundwater lateral fluxes and aquifer balance: \[\begin{equation} \Delta{S_{base,i}} = S_{base, in} - S_{base, out} = \sum_{j}{q_{base, ji}} - \sum_{j}{q_{base, ij}} \end{equation}\]
19.6.2 Overland runoff
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 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. For a given cell \(i\), runoff is defined as the sum of local balance runoff \(Ru_i\) (which includes infiltration excess and saturation excess) and aquifer exfiltration (\(Exf_i\)). This runoff 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 + Exf_i) \cdot q_{ij} \end{equation}\]
Note that a given cell \(j\) can receive overland flux from more than one neighbor.
19.6.3 Watershed surface 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 contributes to watershed runoff exported (\(Exp\)), which is assumed to reach river channels, dams or other surface water bodies.
19.6.4 Losses towards a deeper aquifer
Daily losses from the aquifer towards a deeper aquifer not connected to watershed outlets are simulated by means of a daily maximum loss rate (which can be specified either via the watershed control parameter called deep_aquifer_loss or using an input sf column with the same name) that is compared with the current aquifer content. The minimum of those values is subtracted from the current aquifer content.