hydromodel_spatial: Implementation of Mahrt and Pan two-layer model of soil...

In ilyamaclean/ecohydrotools: hydrogical modelling for ecological applications with R

Description

hydromodel_time is a spatial implementation of the Mahrt and Pan two-layer model of soil hydrology for a point location.

Usage

hydromodel_spatial(dem, rain, cre, sp, basins = NA, theta1 = 0.35,
  theta2 = 0.35, surface = 0, p1 = 0.2, p2 = 0.1, p3 = 0.25,
  z1 = 5, z2 = 95, cn = 82, cover = 0.8, topratio = 0.5,
  timestep = 86400, n2 = 1.1, Ksat2 = 0.5, Trace = TRUE)

Arguments

dem	a raster or matrix of digital elevation data
rain	a vector or 3-dimensional array ofrainfall (mm per time step)
cre	a vector or 3-dimensional array of crop reference evapotranspirations (mm per timestep) as returned by cre.daily() or cre.hourly().
sp	van Genuchten soil paramaters as returned by getsoilparams()
basins	optional raster or matrix of hydrological basins as returned by basindelin. Calculated if not supplied.
theta1	initial soil water fraction of top layer
theta2	initial soil water fraction of bottom layer
surface	initial surface water depth (mm)
p1	a single numeric value of the power adjustment to apply to topographic wetness index values for surface layer (see topdist())
p2	a single numeric value of the power adjustment to apply to topographic wetness index values for sub-surface layer (see topdist())
p3	a single numeric value of the power adjustment to apply to topographic wetness index values for surface water (see topdistw())
z1	a single numeric value, raster or matrix of assumed depth(s) of the top soil layer 1 (cm)
z2	a single numeric value, raster or matrix of assumed depth(s) of bottom soil layer 2 (cm)
cn	a single numeric value, raster or matrix of runoff curve number(s) (see details)
cover	a single numeric value, raster or matrix of fractional of vegetation cover
topratio	a single numeric value, raster or matrix of ratio(s) of root water uptake from top layer relative to bottom layer
timestep	number of seconds in each time step of model run (default one day)
n2	optional single numeric value, raster or matrix of pore size distribution parameter(s) for controlling ground-water seepage (see details)
Ksat2	optional single numeric value, raster or matrix of saturated hydraulic conductivity parameter(s) for controlling ground-water seepage (see details)
Trace	logical indicating whether to produce plots tracking progress

Details

Modified spatial implementation of Mahrt and Pan (1984) two-layer model of soil hydrology. FRactional Soil water content is computed in a thin upper layer for use in calculation of surface evaporation. Water storage is computed for an underlying deeper layer. Precipitation enters the top soil layer, but any precipitation that cannot infiltrate or re-evaporate is specified to be runoff. The rate and direction of exchange of water between the soil layers is determined by soil moisture-dependent hydraulic diffusivity and conductivity and by the difference in soil moisture between the two layers. Evapotranspiration is assumed to occur from the vegetated portion of the basin, and can be partioned between both layers.

In this variant of the model, the spatial distributions of soil and surface water are explicitely accounted for. Hydrological basins are delineated from the dem using basindelin() and the basin volumes, pour points and basins to which overflow would accrue are calculated. Starting with basin with the highest elevation pour point, the model is run iteratively for each basin. Within each time-step soil and surface water are distributed across the basin by topographic wetness index values calculated using topidx(). As such, run-off, hydraulic conductivity and diffusivity are permitted to vary within the basin. Surplus surface water remains within the basin unless the basin volume is exceeded, in which case it is accrued to adjoining basin at the pour point.

Soil moisture-dependent hydraulic diffusivity and conductivity are calculated using the van Genuchten equation (van Genuchten 1980 Soil Science Society of America Journal 44:892-898), and the parameters n, alpha, Ksat, Smax and Smin, obtained using getsoilparams() control the relationship with soil moisture. The paramaters n2 and Ksat2 are equivelent to n and Ksat, but control ground water seepage rates. Setting Ksat2 to zero assumes ground water seepage. Paramaters for the van Genuchten equation should be in units of centimeters and days. Correction for shorter time-scales is automatically applied.

For simplicity, the diffusion of soil water across the bottom of the lower layer is neglected since gradients and resulting fluxes at this depth are generally unimportant, except over longer time-scales. Exceptions include a high water table and the advance of a deep ‘wetting front’. In the latter case, the thickness of the model should be increased.

Surface runoff is dependent on the porosity of the soil and is controlled by a run-off curve number, itself dependent on soil type, land cover and soil condition. Runoff curve numbers can be obtained from here: https://www.wcc.nrcs.usda.gov/ftpref/wntsc/H&H/training/runoff-curve-numbers1.pdf. If add.surplus is set to false (the default), runoff is assumed to leave the basin. If set to true, it is added to rainfall at each time step.

The paramaters rain and cre can be provided as vectors or 3-D arrays. If provides as vectors, spatially uniform rainfall and evapotranspiration is assumed, and number of time-steps over which the model is run is determined the length of these vectors. If provided as 3D arrays, the the first and second dimensions indicated spatial variability, and the third the number of timesteps over which the model is run. The remaining paramaters, with the exception timestep, and Trace can all be provided as single numeric values (no spatial vaiability) or as matrices or rasters (spatial variability).

The model outputs three arrays, each covering the same extent as dem and the third the surface and soil water at each time-step. Note however, that the model relies on the computation of slope angles using the raster package, which sets edge cells to NA. In consequence, the edge cells of the model output are also NA. For the model to work dem must have a planer projection such that x and y values are equidistant.

Value

a list of three arrays: sm1 soil water fractions in the top soil layer, sm2 soil water fractions in the bottom soil layer, surface surface water depth (mm)

See Also

getsoilparams() Dtheta() Ktheta() runoff() soilevap() topdist() topdistw() hydromodel_time()

Examples

# Rainfall and evapotranspiration time series
data(rainevap)
rain <- rainevap$rain
cre <- rainevap$evap
# Get van Genuchten soil paramaters
sp <- getsoilparams("Clay loam")
# ================================================================== #
# Run model for Lizard Peninsula over 365 days assuming clay-loam soils
# and with default settings. Trace is `TRUE` so progress is tracked
# MOdel takes ~ 5 mins to run
# ================================================================== #
hm <- hydromodel_spatial(dtm100m, rain, cre, sp)
# ============== Calculate and plot means  ==================== #
soil1 <- apply(hm$sm1, c(1, 2), mean) # soil moisture in layer 1
soil2 <- apply(hm$sm2, c(1, 2), mean) # soil moisture in layer 2
sfcew <- apply(hm$surface, c(1, 2), mean) # standing water (mm)
plot(if_raster(soil1, dtm100m), main = "Layer 1")
plot(if_raster(soil2, dtm100m), main = "Layer 2")

ilyamaclean/ecohydrotools documentation built on June 10, 2019, 5:45 a.m.